1. All combinations
1.1 Figure3
# Ethanol
#_________________________
rm(list=ls())
#_______________________
# Libraries
#_______________________
library(nortest)
library(writexl)
library(tidyverse)
library(openxlsx)
library(dplyr)
library(Matrix)
library(lme4)
library(ggplot2)
library(car)
library(extrafont)
library(stringr)
library(gridExtra)
library(combinat)
library(tidyr)
library(cowplot)
library(ggpubr)
library(rlang)
library(broom)
library(purrr)
library(magick)
library(ggpmisc) # For adding equation and R² to plot
# Import fonts
# font_import(prompt = FALSE)
loadfonts(device = "win")
get_sharpness_time <- function(data) {
# Loop through the data until the cap value increases significantly
for (i in 2:nrow(data)) {
# Check if there's a notable increase in the Cap value after a steady range
if (data$Cap[i] > 160) {
return(data$Time[i])
}
}
# Return NA if no sharpness time is found
return(NA)
}
calculate_max_time_for_run <- function(data, df2) {
# Extract the first selected run from df2
selected_run <- str_split(df2$combination[1], ", ")[[1]][1]
# Process the data to calculate cumulative time for the selected run
selected_data <- data %>%
mutate(Condition = str_replace(Condition, "-run", "_run")) %>%
filter(Condition == paste0(df2$condition[1], "_", selected_run)) %>%
group_by(Condition) %>%
mutate(Cumulative_Time = cumsum(Time)) %>%
ungroup()
# Return the maximum cumulative time for the selected run
return(max(selected_data$Cumulative_Time))
}
process_selected_data <- function(data, temp_label, df, df2, max_time_templabel) {
# Group and calculate cumulative time
data <- data %>%
group_by(Condition) %>%
mutate(Cumulative_Time = cumsum(Time)) %>%
ungroup() %>%
mutate(Time = Cumulative_Time) %>%
select(-Cumulative_Time)
# Replace "-run" with "_run" in the Condition column to ensure uniformity
data <- data %>%
mutate(Condition = str_replace(Condition, "-run", "_run"))
# Extract Run_Num as a numeric value from Condition
data <- data %>%
mutate(Run_Num = as.numeric(str_extract(Condition, "\\d+$")))
# Ensure df2 and df have consistent column names
df$condition <- factor(df$condition)
# Initialize a list to store plots
plot_list <- list()
# Define condition mapping for Pure Water (0) and Pure Ethanol (100)
condition_map <- c("Pure Water" = 0, "Pure Ethanol" = 100)
for (i in 1:99) {
condition_map[paste0(i, "% Ethanol")] <- i
}
# Define gradient colors for R1, R2, R3
gradient_colors_R1 <- colorRampPalette(c("#6CA6CD", "#000080"))(101) # Light blue to navy blue
gradient_colors_R2 <- colorRampPalette(c("#E396A5", "#8B0000"))(101) # Light red to dark red
gradient_colors_R3 <- colorRampPalette(c("#78D978", "#006400"))(101) # Light green to dark green
# Define fixed colors for R1, R2, and R3 legend (e.g., midpoints or other meaningful colors)
fixed_colors <- c(
"R1" = gradient_colors_R1[51], # Midpoint of R1 gradient (50% Ethanol)
"R2" = gradient_colors_R2[51], # Midpoint of R2 gradient (50% Ethanol)
"R3" = gradient_colors_R3[51] # Midpoint of R3 gradient (50% Ethanol)
)
# Loop over each condition to align data and create plots
for (condition in unique(df2$condition)) {
# Get the corresponding rows for the current condition from df2
df2_subset <- df2 %>% filter(condition == !!condition)
# Initialize a list to store the aligned data for each run
condition_aligned_data <- list()
selected_runs <- unlist(strsplit(as.character(df2_subset$combination[1]), ", "))
# Get the condition percentage from the mapping
condition_percentage <- condition_map[condition]
# Assign colors and linetypes based on the run order
run_linetypes <- c("solid", "dotted", "dotted")
# Iterate through the selected runs
for (i in seq_along(selected_runs)) {
run <- selected_runs[i]
filtered_data <- data %>% filter(Condition == paste(condition, run, sep = "_"))
# Get the sharpness time and align the data
sharpness_time <- get_sharpness_time(filtered_data)
data_run <- filtered_data %>%
mutate(Time_aligned = Time - sharpness_time) %>%
mutate(Time_aligned = pmax(Time_aligned, 0)) %>% # Set negative values to 0
mutate(Cap_aligned = Cap) # Retain Cap values
# Assign gradient color based on the condition percentage
gradient_color <- switch(i,
gradient_colors_R1[condition_percentage + 1],
gradient_colors_R2[condition_percentage + 1],
gradient_colors_R3[condition_percentage + 1]
)
condition_aligned_data[[i]] <- data_run %>%
mutate(Run = paste0("R", i), Run_Color = gradient_color, Run_Linetype = run_linetypes[i])
}
# Combine all aligned data for this condition
aligned_data_combined <- bind_rows(condition_aligned_data)
# Filter data based on max time
aligned_data_combined <- aligned_data_combined %>%
filter(Time_aligned <= max_time_templabel)
# Create the plot
plot <- aligned_data_combined %>%
ggplot(aes(x = Time_aligned, y = Cap_aligned, group = Run, color = Run_Color, linetype = Run_Linetype)) +
geom_line(size = 1.2) +
scale_color_identity(
guide = guide_legend(override.aes = list(color = fixed_colors, linetype = run_linetypes)),
labels = names(fixed_colors)
) +
scale_linetype_identity() +
labs(
x = "Time(s)",
y = "C(fF)",
color = "Runs",
linetype = "Runs"
) +
theme_minimal() +
theme(
axis.text = element_text(size = 20, family = "Times New Roman", face = "bold"),
axis.text.x = element_text(margin = margin(t = 10)),
axis.title.x = element_text(size = 20, family = "Times New Roman", face = "bold", margin = margin(t = 0, r = 40, b = 0, l = 10)),
axis.title.y = element_text(size = 20, family = "Times New Roman", face = "bold", margin = margin(t = 0, r = 40, b = 0, l = 10)),
legend.text = element_text(size = 15, family = "Times New Roman", face = "bold"),
legend.title = element_text(size = 15, family = "Times New Roman", face = "bold"),
legend.position = c(0.7, 0.5),
legend.direction = "vertical",
legend.spacing.y = unit(0.002, 'cm'),
legend.key.size = unit(0.8, 'lines'),
panel.background = element_blank(),
plot.background = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.border = element_rect(color = "black", fill = NA)
) +
scale_x_continuous(
limits = c(0, 850),
breaks = seq(0, 800, by = 400)) +
coord_cartesian(ylim = c(130, 180)) +
scale_y_continuous(breaks = c(130, 150, 170))
plot_list[[condition]] <- plot
}
# Combine the plots into a grid
combined_plot <- cowplot::plot_grid(
plotlist = lapply(plot_list, function(p) {
p
}),
nrow = 11,
ncol = 1,
align = "v"
)
# Save the combined plot
ggsave(
filename = paste0("Combined_Plot_E", temp_label, ".png"),
plot = combined_plot,
height = 20,
width = 2.5
)
}
data_25C <- read.xlsx("reshape_Water-Ethanol.xlsx", sheet = "25C", colNames = TRUE)
data_40C <- read.xlsx("reshape_Water-Ethanol.xlsx", sheet = "40C", colNames = TRUE)
data_50C <- read.xlsx("reshape_Water-Ethanol.xlsx", sheet = "50C", colNames = TRUE)
data_60C <- read.xlsx("reshape_Water-Ethanol.xlsx", sheet = "60C", colNames = TRUE)
df_25C <- read.xlsx("timeline_results__ethanol_water25C.xlsx", colNames = TRUE)
df_40C <- read.xlsx("timeline_results__ethanol_water40C.xlsx", colNames = TRUE)
df2_40C <- read.xlsx("mean_sd_combinations_and_selected_runs_ethanol_water40C.xlsx", sheet = "Best_Combinations", colNames = TRUE)
df2_25C <- read.xlsx("mean_sd_combinations_and_selected_runs_ethanol_water25C.xlsx", sheet = "Best_Combinations", colNames = TRUE)
df_50C <- read.xlsx("timeline_results__ethanol_water50C.xlsx", colNames = TRUE)
df2_50C <- read.xlsx("mean_sd_combinations_and_selected_runs_ethanol_water50C.xlsx", sheet = "Best_Combinations", colNames = TRUE)
df_60C <- read.xlsx("timeline_results__ethanol_water60C.xlsx", colNames = TRUE)
df2_60C <- read.xlsx("mean_sd_combinations_and_selected_runs_ethanol_water60C.xlsx", sheet = "Best_Combinations", colNames = TRUE)
# Calculate the maximum time for each temperature's selected run
max_time_25C <- calculate_max_time_for_run(data_25C, df2_25C)
max_time_40C <- calculate_max_time_for_run(data_40C, df2_40C)
max_time_50C <- calculate_max_time_for_run(data_50C, df2_50C)
max_time_60C <- calculate_max_time_for_run(data_60C, df2_60C)
# Apply the process_data function to each temperature
process_selected_data(data_25C, "25C", df_25C, df2_25C, max_time_25C)
process_selected_data(data_40C, "40C", df_40C, df2_40C, max_time_40C)
process_selected_data(data_50C, "50C", df_50C, df2_50C, max_time_50C)
process_selected_data(data_60C, "60C", df_60C, df2_60C, max_time_60C)
############selected curves#############
############selected curves#############
# Function to calculate the maximum cumulative time for the selected run of a temperature
calculate_max_time_for_run <- function(data, df2) {
# Extract the first selected run from df2
selected_run <- str_split(df2$combination[1], ", ")[[1]][1]
# Process the data to calculate cumulative time for the selected run
selected_data <- data %>%
mutate(Condition = str_replace(Condition, "-run", "_run")) %>%
filter(Condition == paste0(df2$condition[1], "_", selected_run)) %>%
group_by(Condition) %>%
mutate(Cumulative_Time = cumsum(Time)) %>%
ungroup()
# Return the maximum cumulative time for the selected run
return(max(selected_data$Cumulative_Time))
}
# Calculate the maximum time for each temperature's selected run
max_time_25C <- calculate_max_time_for_run(data_25C, df2_25C)
max_time_40C <- calculate_max_time_for_run(data_40C, df2_40C)
max_time_50C <- calculate_max_time_for_run(data_50C, df2_50C)
max_time_60C <- calculate_max_time_for_run(data_60C, df2_60C)
# Function to process and visualize data for each temperature using the specific max time
process_selected_data <- function(data, temp_label, df, df2, max_time_templabel) {
data <- data %>%
group_by(Condition) %>%
mutate(Cumulative_Time = cumsum(Time)) %>%
ungroup() %>%
mutate(Time = Cumulative_Time) %>%
select(-Cumulative_Time) %>%
mutate(Condition = str_replace(Condition, "-run", "_run"))
# Ensure 'df' is structured properly
df$condition <- factor(df$condition)
# Preparing to store selected data
selected_data <- data.frame()
# Create ethanol labels for each condition depending on the temperature
if (temp_label == "25C") {
ethanol_labels <- c(
"Pure Water_run10" = "0",
"10% Ethanol_run5" = "10",
"20% Ethanol_run4" = "20",
"30% Ethanol_run5" = "30",
"40% Ethanol_run5" = "40",
"50% Ethanol_run1" = "50",
"60% Ethanol_run2" = "60",
"70% Ethanol_run2" = "70",
"80% Ethanol_run1" = "80",
"90% Ethanol_run3" = "90",
"Pure Ethanol_run3" = "100"
)
} else {
ethanol_labels <- c(
"Pure Water_run1" = "0",
"10% Ethanol_run1" = "10",
"20% Ethanol_run1" = "20",
"30% Ethanol_run1" = "30",
"40% Ethanol_run1" = "40",
"50% Ethanol_run1" = "50",
"60% Ethanol_run1" = "60",
"70% Ethanol_run1" = "70",
"80% Ethanol_run1" = "80",
"90% Ethanol_run1" = "90",
"Pure Ethanol_run1" = "100"
)
}
# Process each condition from df2
unique_conditions <- unique(df2$condition)
for (condition in unique_conditions) {
df2_subset <- df2 %>% filter(condition == !!condition)
selected_runs <- unlist(strsplit(as.character(df2_subset$combination[1]), ", "))
run <- selected_runs[1] # Ensure 'run' exists in our data columns
filtered_data <- data %>%
filter(Condition == paste(condition, run, sep = "_"))
# Get the sharpness time for aligning data
sharpness_time <- get_sharpness_time(filtered_data)
data_run <- filtered_data %>%
mutate(Time_aligned = pmax(Time - sharpness_time, 0)) %>%
mutate(Cap_aligned = Cap)
selected_data <- bind_rows(selected_data, data_run %>% select(Time_aligned, Cap_aligned, Condition))
}
# Determine a common starting value, e.g., the median of the first capacitance values of each condition
common_start_value <- median(selected_data$Cap_aligned[selected_data$Time_aligned == 0])
# Function to adjust the capacitance values so that all start from the common start value
adjust_capacitance <- function(data, common_value) {
data %>%
group_by(Condition) %>%
mutate(adjustment = common_value - first(Cap_aligned)) %>%
ungroup() %>%
mutate(Cap_aligned = Cap_aligned + adjustment)
}
# Apply the adjustment function
adjusted_data <- adjust_capacitance(selected_data, common_start_value)
# Truncate data at the max time for the current temperature
adjusted_data <- adjusted_data %>%
filter(Time_aligned <= max_time_templabel)
# Update the Condition column to reflect ethanol percentages
adjusted_data <- adjusted_data %>%
mutate(Condition = ethanol_labels[Condition])
adjusted_data$Condition <- factor(adjusted_data$Condition, levels = c("0", "10", "20", "30", "40", "50", "60", "70", "80", "90", "100"))
# Generate gradient colors from light blue to navy blue for 11 conditions (0% to 100%)
color_gradient <- colorRampPalette(c("#6CA6CD", "navyblue"))(11)
# Plot the adjusted data with custom color and legend settings
selected_plot <- ggplot(adjusted_data, aes(x = Time_aligned, y = Cap_aligned, color = Condition)) +
geom_line(size = 1) +
labs(
x = "Time(s)",
y = "C(fF)",
color = paste("E_W(%),", temp_label)) + # Adding temp_label to the legend title
theme_minimal() +
theme(
axis.text = element_text(size = 20, family = "Times New Roman", face = "bold"),
axis.text.x = element_text(margin = margin(t = 10)), # Add margin to x-axis text
axis.title.x = element_text(size = 20, family = "Times New Roman", face = "bold", margin = margin(t = 0, r = 10, b = 0, l = 5)),
axis.title.y = element_text(size = 20, family = "Times New Roman", face = "bold", margin = margin(t = 0, r = 10, b = 0, l = 5)),
legend.text = element_text(size = 12, family = "Times New Roman", face = "bold"),
legend.title = element_text(size = 12, family = "Times New Roman", face = "bold"),
legend.position = "none",
legend.direction = "vertical",
legend.spacing.y = unit(0.002, 'cm'),
legend.key.size = unit(0.8, 'lines'),
panel.background = element_blank(),
plot.background = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.border = element_rect(color = "black", fill = NA)
) +
# Apply scale_color_manual with dynamically generated gradient colors
scale_color_manual(
values = color_gradient
)+
scale_x_continuous(
limits = c(0, 850), # Slightly extend the limit
breaks = seq(0, 800, by = 400)) +
coord_cartesian(ylim = c(130, 180)) +
scale_y_continuous(breaks = c(130, 150, 170))
# Save the plot
ggsave(filename = paste0("Selected_Curves_E", temp_label, ".png"), plot = selected_plot, height = 2, width = 2.5)
}
# Apply the process_selected_data function to the data of each temperature using its specific max time
process_selected_data(data_25C, "25C", df_25C, df2_25C, max_time_25C)
process_selected_data(data_40C, "40C", df_40C, df2_40C, max_time_40C)
process_selected_data(data_50C, "50C", df_50C, df2_50C, max_time_50C)
process_selected_data(data_60C, "60C", df_60C, df2_60C, max_time_60C)
############
combine_and_align <- function(upper_image_path, lower_image_path, output_path) {
# Load images
upper_image <- image_read(upper_image_path)
lower_image <- image_read(lower_image_path)
# Trim the images to calculate the content width
trimmed_upper_image <- image_trim(upper_image)
trimmed_lower_image <- image_trim(lower_image)
# Get the width of the trimmed images
trimmed_upper_width <- image_info(trimmed_upper_image)$width
trimmed_lower_width <- image_info(trimmed_lower_image)$width
# Calculate the left padding (difference between original width and trimmed width)
left_padding_upper <- image_info(upper_image)$width - trimmed_upper_width
left_padding_lower <- image_info(lower_image)$width - trimmed_lower_width
# Determine the maximum width for alignment
aligned_width <- max(image_info(upper_image)$width, image_info(lower_image)$width)
# Adjust both images with consistent left padding and total width
aligned_upper_image <- image_extent(upper_image, paste0(aligned_width, "x"), gravity = "west")
aligned_lower_image <- image_extent(lower_image, paste0(aligned_width, "x"), gravity = "west")
# Stack images vertically
combined_image <- image_append(c(aligned_upper_image, aligned_lower_image), stack = TRUE)
# Save the final aligned image
image_write(combined_image, path = output_path)
}
# Use the function with your specific files
combine_and_align("Combined_Plot_E25C.png", "Selected_Curves_E25C.png", "Combined_Plot_E25C.png")
combine_and_align("Combined_Plot_E40C.png", "Selected_Curves_E40C.png", "Combined_Plot_E40C.png")
combine_and_align("Combined_Plot_E50C.png", "Selected_Curves_E50C.png", "Combined_Plot_E50C.png")
combine_and_align("Combined_Plot_E60C.png", "Selected_Curves_E60C.png", "Combined_Plot_E60C.png")
#__________________________
# Methanol
#_________________________
rm(list=ls())
#_______________________
# Libraries
#_______________________
library(nortest)
library(writexl)
library(tidyverse)
library(openxlsx)
library(dplyr)
library(Matrix)
library(lme4)
library(ggplot2)
library(car)
library(extrafont)
library(stringr)
library(gridExtra)
library(combinat)
library(tidyr)
library(cowplot)
library(ggpubr)
library(rlang)
library(broom)
library(purrr)
library(magick)
library(ggpmisc) # For adding equation and R² to plot
# Import fonts
# font_import(prompt = FALSE)
loadfonts(device = "win")
get_sharpness_time <- function(data) {
# Loop through the data until the cap value increases significantly
for (i in 2:nrow(data)) {
# Check if there's a notable increase in the Cap value after a steady range
if (data$Cap[i] > 160) {
return(data$Time[i])
}
}
# Return NA if no sharpness time is found
return(NA)
}
calculate_max_time_for_run <- function(data, df2) {
# Extract the first selected run from df2
selected_run <- str_split(df2$combination[1], ", ")[[1]][1]
# Process the data to calculate cumulative time for the selected run
selected_data <- data %>%
mutate(Condition = str_replace(Condition, "-run", "_run")) %>%
filter(Condition == paste0(df2$condition[1], "_", selected_run)) %>%
group_by(Condition) %>%
mutate(Cumulative_Time = cumsum(Time)) %>%
ungroup()
# Return the maximum cumulative time for the selected run
return(max(selected_data$Cumulative_Time))
}
process_selected_data <- function(data, temp_label, df, df2, max_time_templabel) {
# Group and calculate cumulative time
data <- data %>%
group_by(Condition) %>%
mutate(Cumulative_Time = cumsum(Time)) %>%
ungroup() %>%
mutate(Time = Cumulative_Time) %>%
select(-Cumulative_Time)
# Replace "-run" with "_run" in the Condition column to ensure uniformity
data <- data %>%
mutate(Condition = str_replace(Condition, "-run", "_run"))
# Extract Run_Num as a numeric value from Condition
data <- data %>%
mutate(Run_Num = as.numeric(str_extract(Condition, "\\d+$")))
# Ensure df2 and df have consistent column names
df$condition <- factor(df$condition)
# Initialize a list to store plots
plot_list <- list()
# Define condition mapping for Pure Water (0) and Pure Methanol (100)
condition_map <- c("Pure Water" = 0, "Pure Methanol" = 100)
for (i in 1:99) {
condition_map[paste0(i, "% Methanol")] <- i
}
# Define gradient colors for R1, R2, R3
gradient_colors_R1 <- colorRampPalette(c("#D2B48C", "#8B4513"))(101) # Light blue to navy blue
gradient_colors_R2 <- colorRampPalette(c("#E396A5", "#8B0000"))(101) # Light red to dark red
gradient_colors_R3 <- colorRampPalette(c("#78D978", "#006400"))(101) # Light green to dark green
# Define fixed colors for R1, R2, and R3 legend (e.g., midpoints or other meaningful colors)
fixed_colors <- c(
"R1" = gradient_colors_R1[51], # Midpoint of R1 gradient (50% Methanol)
"R2" = gradient_colors_R2[51], # Midpoint of R2 gradient (50% Methanol)
"R3" = gradient_colors_R3[51] # Midpoint of R3 gradient (50% Methanol)
)
# Loop over each condition to align data and create plots
for (condition in unique(df2$condition)) {
# Get the corresponding rows for the current condition from df2
df2_subset <- df2 %>% filter(condition == !!condition)
# Initialize a list to store the aligned data for each run
condition_aligned_data <- list()
selected_runs <- unlist(strsplit(as.character(df2_subset$combination[1]), ", "))
# Get the condition percentage from the mapping
condition_percentage <- condition_map[condition]
# Assign colors and linetypes based on the run order
run_linetypes <- c("solid", "dotted", "dotted")
# Iterate through the selected runs
for (i in seq_along(selected_runs)) {
run <- selected_runs[i]
filtered_data <- data %>% filter(Condition == paste(condition, run, sep = "_"))
# Get the sharpness time and align the data
sharpness_time <- get_sharpness_time(filtered_data)
data_run <- filtered_data %>%
mutate(Time_aligned = Time - sharpness_time) %>%
mutate(Time_aligned = pmax(Time_aligned, 0)) %>% # Set negative values to 0
mutate(Cap_aligned = Cap) # Retain Cap values
# Assign gradient color based on the condition percentage
gradient_color <- switch(i,
gradient_colors_R1[condition_percentage + 1],
gradient_colors_R2[condition_percentage + 1],
gradient_colors_R3[condition_percentage + 1]
)
condition_aligned_data[[i]] <- data_run %>%
mutate(Run = paste0("R", i), Run_Color = gradient_color, Run_Linetype = run_linetypes[i])
}
# Combine all aligned data for this condition
aligned_data_combined <- bind_rows(condition_aligned_data)
# Filter data based on max time
aligned_data_combined <- aligned_data_combined %>%
filter(Time_aligned <= max_time_templabel)
# Create the plot
plot <- aligned_data_combined %>%
ggplot(aes(x = Time_aligned, y = Cap_aligned, group = Run, color = Run_Color, linetype = Run_Linetype)) +
geom_line(size = 1.2) +
scale_color_identity(
guide = guide_legend(override.aes = list(color = fixed_colors, linetype = run_linetypes)),
labels = names(fixed_colors)
) +
scale_linetype_identity() +
labs(
x = "Time(s)",
y = "C(fF)",
color = "Runs",
linetype = "Runs"
) +
theme_minimal() +
theme(
axis.text = element_text(size = 20, family = "Times New Roman", face = "bold"),
axis.text.x = element_text(margin = margin(t = 10)),
axis.title.x = element_text(size = 20, family = "Times New Roman", face = "bold", margin = margin(t = 0, r = 40, b = 0, l = 10)),
axis.title.y = element_text(size = 20, family = "Times New Roman", face = "bold", margin = margin(t = 0, r = 40, b = 0, l = 10)),
legend.text = element_text(size = 15, family = "Times New Roman", face = "bold"),
legend.title = element_text(size = 15, family = "Times New Roman", face = "bold"),
legend.position = c(0.7, 0.5),
legend.direction = "vertical",
legend.spacing.y = unit(0.002, 'cm'),
legend.key.size = unit(0.8, 'lines'),
panel.background = element_blank(),
plot.background = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.border = element_rect(color = "black", fill = NA)
) +
scale_x_continuous(
limits = c(0, 850),
breaks = seq(0, 800, by = 400)) +
coord_cartesian(ylim = c(130, 180)) +
scale_y_continuous(breaks = c(130, 150, 170))
plot_list[[condition]] <- plot
}
# Combine the plots into a grid
combined_plot <- cowplot::plot_grid(
plotlist = lapply(plot_list, function(p) {
p
}),
nrow = 11,
ncol = 1,
align = "v"
)
# Save the combined plot
ggsave(
filename = paste0("Combined_Plot_M", temp_label, ".png"),
plot = combined_plot,
height = 20,
width = 2.5
)
}
data_25C <- read.xlsx("reshape_Water-Methanol.xlsx", sheet = "25C", colNames = TRUE)
data_40C <- read.xlsx("reshape_Water-Methanol.xlsx", sheet = "40C", colNames = TRUE)
data_50C <- read.xlsx("reshape_Water-Methanol.xlsx", sheet = "50C", colNames = TRUE)
data_60C <- read.xlsx("reshape_Water-Methanol.xlsx", sheet = "60C", colNames = TRUE)
df_25C <- read.xlsx("timeline_results__Methanol_water25C.xlsx", colNames = TRUE)
df_40C <- read.xlsx("timeline_results__Methanol_water40C.xlsx", colNames = TRUE)
df2_40C <- read.xlsx("mean_sd_combinations_and_selected_runs_Methanol_water40C.xlsx", sheet = "Best_Combinations", colNames = TRUE)
df2_25C <- read.xlsx("mean_sd_combinations_and_selected_runs_Methanol_water25C.xlsx", sheet = "Best_Combinations", colNames = TRUE)
df_50C <- read.xlsx("timeline_results__Methanol_water50C.xlsx", colNames = TRUE)
df2_50C <- read.xlsx("mean_sd_combinations_and_selected_runs_Methanol_water50C.xlsx", sheet = "Best_Combinations", colNames = TRUE)
df_60C <- read.xlsx("timeline_results__Methanol_water60C.xlsx", colNames = TRUE)
df2_60C <- read.xlsx("mean_sd_combinations_and_selected_runs_Methanol_water60C.xlsx", sheet = "Best_Combinations", colNames = TRUE)
# Calculate the maximum time for each temperature's selected run
max_time_25C <- calculate_max_time_for_run(data_25C, df2_25C)
max_time_40C <- calculate_max_time_for_run(data_40C, df2_40C)
max_time_50C <- calculate_max_time_for_run(data_50C, df2_50C)
max_time_60C <- calculate_max_time_for_run(data_60C, df2_60C)
# Apply the process_data function to each temperature
process_selected_data(data_25C, "25C", df_25C, df2_25C, max_time_25C)
process_selected_data(data_40C, "40C", df_40C, df2_40C, max_time_40C)
process_selected_data(data_50C, "50C", df_50C, df2_50C, max_time_50C)
process_selected_data(data_60C, "60C", df_60C, df2_60C, max_time_60C)
############selected curves#############
# Function to calculate the maximum cumulative time for the selected run of a temperature
calculate_max_time_for_run <- function(data, df2) {
# Extract the first selected run from df2
selected_run <- str_split(df2$combination[1], ", ")[[1]][1]
# Process the data to calculate cumulative time for the selected run
selected_data <- data %>%
mutate(Condition = str_replace(Condition, "-run", "_run")) %>%
filter(Condition == paste0(df2$condition[1], "_", selected_run)) %>%
group_by(Condition) %>%
mutate(Cumulative_Time = cumsum(Time)) %>%
ungroup()
# Return the maximum cumulative time for the selected run
return(max(selected_data$Cumulative_Time))
}
# Calculate the maximum time for each temperature's selected run
max_time_25C <- calculate_max_time_for_run(data_25C, df2_25C)
max_time_40C <- calculate_max_time_for_run(data_40C, df2_40C)
max_time_50C <- calculate_max_time_for_run(data_50C, df2_50C)
max_time_60C <- calculate_max_time_for_run(data_60C, df2_60C)
# Function to process and visualize data for each temperature using the specific max time
process_selected_data <- function(data, temp_label, df, df2, max_time_templabel) {
data <- data %>%
group_by(Condition) %>%
mutate(Cumulative_Time = cumsum(Time)) %>%
ungroup() %>%
mutate(Time = Cumulative_Time) %>%
select(-Cumulative_Time) %>%
mutate(Condition = str_replace(Condition, "-run", "_run"))
# Ensure 'df' is structured properly
df$condition <- factor(df$condition)
# Preparing to store selected data
selected_data <- data.frame()
if (temp_label == "25C") {
Methanol_labels <- c(
"Pure Water_run7" = "0",
"10% Methanol_run3" = "10",
"20% Methanol_run2" = "20",
"30% Methanol_run1" = "30",
"40% Methanol_run2" = "40",
"50% Methanol_run1" = "50",
"60% Methanol_run3" = "60",
"70% Methanol_run4" = "70",
"80% Methanol_run4" = "80",
"90% Methanol_run1" = "90",
"Pure Methanol_run1" = "100"
)
} else if (temp_label == "40C") {
Methanol_labels <- c(
"Pure Water_run2" = "0",
"10% Methanol_run1" = "10",
"20% Methanol_run1" = "20",
"30% Methanol_run1" = "30",
"40% Methanol_run1" = "40",
"50% Methanol_run1" = "50",
"60% Methanol_run1" = "60",
"70% Methanol_run1" = "70",
"80% Methanol_run1" = "80",
"90% Methanol_run1" = "90",
"Pure Methanol_run1" = "100"
)
} else if (temp_label == "60C") {
Methanol_labels <- c(
"Pure Water_run4" = "0",
"10% Methanol_run1" = "10",
"20% Methanol_run1" = "20",
"30% Methanol_run1" = "30",
"40% Methanol_run1" = "40",
"50% Methanol_run1" = "50",
"60% Methanol_run1" = "60",
"70% Methanol_run1" = "70",
"80% Methanol_run1" = "80",
"90% Methanol_run1" = "90",
"Pure Methanol_run1" = "100"
)
} else {
Methanol_labels <- c(
"Pure Water_run1" = "0",
"10% Methanol_run1" = "10",
"20% Methanol_run1" = "20",
"30% Methanol_run1" = "30",
"40% Methanol_run1" = "40",
"50% Methanol_run1" = "50",
"60% Methanol_run1" = "60",
"70% Methanol_run1" = "70",
"80% Methanol_run1" = "80",
"90% Methanol_run1" = "90",
"Pure Methanol_run1" = "100"
)
}
# Process each condition from df2
unique_conditions <- unique(df2$condition)
for (condition in unique_conditions) {
df2_subset <- df2 %>% filter(condition == !!condition)
selected_runs <- unlist(strsplit(as.character(df2_subset$combination[1]), ", "))
run <- selected_runs[1] # Ensure 'run' exists in our data columns
filtered_data <- data %>%
filter(Condition == paste(condition, run, sep = "_"))
filtered_data <- filtered_data %>%
filter(!is.na(Cap))
# Get the sharpness time for aligning data
sharpness_time <- get_sharpness_time(filtered_data)
data_run <- filtered_data %>%
mutate(Time_aligned = pmax(Time - sharpness_time, 0)) %>%
mutate(Cap_aligned = Cap)
selected_data <- bind_rows(selected_data, data_run %>% select(Time_aligned, Cap_aligned, Condition))
}
# Determine a common starting value, e.g., the median of the first capacitance values of each condition
common_start_value <- median(selected_data$Cap_aligned[selected_data$Time_aligned == 0])
# Function to adjust the capacitance values so that all start from the common start value
adjust_capacitance <- function(data, common_value) {
data %>%
group_by(Condition) %>%
mutate(adjustment = common_value - first(Cap_aligned)) %>%
ungroup() %>%
mutate(Cap_aligned = Cap_aligned + adjustment)
}
# Apply the adjustment function
adjusted_data <- adjust_capacitance(selected_data, common_start_value)
# Truncate data at the max time for the current temperature
adjusted_data <- adjusted_data %>%
filter(Time_aligned <= max_time_templabel)
# Update the Condition column to reflect ethanol percentages
adjusted_data <- adjusted_data %>%
mutate(Condition = Methanol_labels[Condition])
adjusted_data$Condition <- factor(adjusted_data$Condition, levels = c("0", "10", "20", "30", "40", "50", "60", "70", "80", "90", "100"))
# Generate gradient colors from light blue to navy blue for 11 conditions (0% to 100%)
color_gradient <- colorRampPalette(c("#D2B48C", "#8B4513"))(11)
# Plot the adjusted data with custom color and legend settings
selected_plot <- ggplot(adjusted_data, aes(x = Time_aligned, y = Cap_aligned, color = Condition)) +
geom_line(size = 1) +
labs(
x = "Time(s)",
y = "C(fF)",
color = paste("M_W(%),", temp_label)) + # Adding temp_label to the legend title
theme_minimal() +
theme(
axis.text = element_text(size = 20, family = "Times New Roman", face = "bold"),
axis.text.x = element_text(margin = margin(t = 10)), # Add margin to x-axis text
axis.title.x = element_text(size = 20, family = "Times New Roman", face = "bold", margin = margin(t = 0, r = 10, b = 0, l = 5)),
axis.title.y = element_text(size = 20, family = "Times New Roman", face = "bold", margin = margin(t = 0, r = 10, b = 0, l = 5)),
legend.text = element_text(size = 12, family = "Times New Roman", face = "bold"),
legend.title = element_text(size = 12, family = "Times New Roman", face = "bold"),
legend.position = "none",
legend.direction = "vertical",
legend.spacing.y = unit(0.002, 'cm'),
legend.key.size = unit(0.8, 'lines'),
panel.background = element_blank(),
plot.background = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.border = element_rect(color = "black", fill = NA)
) +
# Apply scale_color_manual with dynamically generated gradient colors
scale_color_manual(
values = color_gradient
)+
scale_x_continuous(
limits = c(0, 850), # Slightly extend the limit
breaks = seq(0, 800, by = 400)) +
coord_cartesian(
ylim = if (temp_label == "60C") c(130, 190) else c(130, 180)
) +
scale_y_continuous(
breaks = if (temp_label == "60C") c(130, 150, 180) else c(130, 150, 170)
)
# Save the plot
ggsave(filename = paste0("Selected_Curves_M", temp_label, ".png"), plot = selected_plot, height = 2, width = 2.5)
}
# Apply the process_selected_data function to the data of each temperature using its specific max time
process_selected_data(data_25C, "25C", df_25C, df2_25C, max_time_25C)
process_selected_data(data_40C, "40C", df_40C, df2_40C, max_time_40C)
process_selected_data(data_50C, "50C", df_50C, df2_50C, max_time_50C)
process_selected_data(data_60C, "60C", df_60C, df2_60C, max_time_60C)
############
combine_and_align <- function(upper_image_path, lower_image_path, output_path) {
# Load images
upper_image <- image_read(upper_image_path)
lower_image <- image_read(lower_image_path)
# Trim the images to calculate the content width
trimmed_upper_image <- image_trim(upper_image)
trimmed_lower_image <- image_trim(lower_image)
# Get the width of the trimmed images
trimmed_upper_width <- image_info(trimmed_upper_image)$width
trimmed_lower_width <- image_info(trimmed_lower_image)$width
# Calculate the left padding (difference between original width and trimmed width)
left_padding_upper <- image_info(upper_image)$width - trimmed_upper_width
left_padding_lower <- image_info(lower_image)$width - trimmed_lower_width
# Determine the maximum width for alignment
aligned_width <- max(image_info(upper_image)$width, image_info(lower_image)$width)
# Adjust both images with consistent left padding and total width
aligned_upper_image <- image_extent(upper_image, paste0(aligned_width, "x"), gravity = "west")
aligned_lower_image <- image_extent(lower_image, paste0(aligned_width, "x"), gravity = "west")
# Stack images vertically
combined_image <- image_append(c(aligned_upper_image, aligned_lower_image), stack = TRUE)
# Save the final aligned image
image_write(combined_image, path = output_path)
}
# Use the function with your specific files
combine_and_align("Combined_Plot_M25C.png", "Selected_Curves_M25C.png", "Combined_Plot_M25C.png")
combine_and_align("Combined_Plot_M40C.png", "Selected_Curves_M40C.png", "Combined_Plot_M40C.png")
combine_and_align("Combined_Plot_M50C.png", "Selected_Curves_M50C.png", "Combined_Plot_M50C.png")
combine_and_align("Combined_Plot_M60C.png", "Selected_Curves_M60C.png", "Combined_Plot_M60C.png")
############__________________________
# Ethanol_Methanol
#_________________________
############################################
#_________________________
rm(list=ls())
#_______________________
# Libraries
#_______________________
library(nortest)
library(writexl)
library(tidyverse)
library(openxlsx)
library(dplyr)
library(Matrix)
library(lme4)
library(ggplot2)
library(car)
library(extrafont)
library(stringr)
library(gridExtra)
library(combinat)
library(tidyr)
library(cowplot)
library(ggpubr)
library(rlang)
library(broom)
library(purrr)
library(magick)
library(ggpmisc) # For adding equation and R² to plot
# Import fonts
# font_import(prompt = FALSE)
loadfonts(device = "win")
get_sharpness_time <- function(data) {
# Loop through the data until the cap value increases significantly
for (i in 2:nrow(data)) {
# Check if there's a notable increase in the Cap value after a steady range
if (data$Cap[i] > 160) {
return(data$Time[i])
}
}
# Return NA if no sharpness time is found
return(NA)
}
data_25C <- read.xlsx("reshape_Methanol-Ethanol.xlsx", sheet = "25C", colNames = TRUE)
df_25C <- read.xlsx("timeline_results__Methanol_Ethanol25C.xlsx", colNames = TRUE)
df2_25C <- read.xlsx("mean_sd_combinations_and_selected_runs_methanol_ethanol25C.xlsx", sheet = "Best_Combinations", colNames = TRUE)
calculate_max_time_for_run <- function(data, df2) {
# Extract the first selected run from df2
selected_run <- str_split(df2$combination[1], ", ")[[1]][1]
# Process the data to calculate cumulative time for the selected run
selected_data <- data %>%
mutate(Condition = str_replace(Condition, "-run", "_run")) %>%
filter(Condition == paste0(df2$condition[1], "_", selected_run)) %>%
group_by(Condition) %>%
mutate(Cumulative_Time = cumsum(Time)) %>%
ungroup()
# Return the maximum cumulative time for the selected run
return(max(selected_data$Cumulative_Time))
}
# Calculate the maximum time for each temperature's selected run
max_time_25C <- calculate_max_time_for_run(data_25C, df2_25C)
process_selected_data <- function(data, temp_label, df, df2, max_time_templabel) {
# Group and calculate cumulative time
data <- data %>%
group_by(Condition) %>%
mutate(Cumulative_Time = cumsum(Time)) %>%
ungroup() %>%
mutate(Time = Cumulative_Time) %>%
select(-Cumulative_Time)
# Replace "-run" with "_run" in the Condition column to ensure uniformity
data <- data %>%
mutate(Condition = str_replace(Condition, "-run", "_run"))
# Extract Run_Num as a numeric value from Condition
data <- data %>%
mutate(Run_Num = as.numeric(str_extract(Condition, "\\d+$")))
# Ensure df2 and df have consistent column names
df$condition <- factor(df$condition)
# Initialize a list to store plots
plot_list <- list()
# Define condition mapping for Pure Water (0) and Pure Methanol (100)
condition_map <- c("Pure Methanol" = 100)
for (i in 0:99) {
condition_map[paste0(i, "% Methanol")] <- i
}
# Define gradient colors for R1, R2, R3
gradient_colors_R1 <- colorRampPalette(c("#66a61e", "#00441b"))(101) # Light blue to navy blue
gradient_colors_R2 <- colorRampPalette(c("#E396A5", "#8B0000"))(101) # Light red to dark red
gradient_colors_R3 <- colorRampPalette(c("#78D978", "#006400"))(101) # Light green to dark green
# Define fixed colors for R1, R2, and R3 legend (e.g., midpoints or other meaningful colors)
fixed_colors <- c(
"R1" = gradient_colors_R1[51], # Midpoint of R1 gradient (50% Methanol)
"R2" = gradient_colors_R2[51], # Midpoint of R2 gradient (50% Methanol)
"R3" = gradient_colors_R3[51] # Midpoint of R3 gradient (50% Methanol)
)
# Loop over each condition to align data and create plots
for (condition in unique(df2$condition)) {
# Get the corresponding rows for the current condition from df2
df2_subset <- df2 %>% filter(condition == !!condition)
# Initialize a list to store the aligned data for each run
condition_aligned_data <- list()
selected_runs <- unlist(strsplit(as.character(df2_subset$combination[1]), ", "))
# Get the condition percentage from the mapping
condition_percentage <- condition_map[condition]
# Assign colors and linetypes based on the run order
run_linetypes <- c("solid", "dotted", "dotted")
# Iterate through the selected runs
for (i in seq_along(selected_runs)) {
run <- selected_runs[i]
filtered_data <- data %>% filter(Condition == paste(condition, run, sep = "_"))
# Get the sharpness time and align the data
sharpness_time <- get_sharpness_time(filtered_data)
data_run <- filtered_data %>%
mutate(Time_aligned = Time - sharpness_time) %>%
mutate(Time_aligned = pmax(Time_aligned, 0)) %>% # Set negative values to 0
mutate(Cap_aligned = Cap) # Retain Cap values
# Assign gradient color based on the condition percentage
gradient_color <- switch(i,
gradient_colors_R1[condition_percentage + 1],
gradient_colors_R2[condition_percentage + 1],
gradient_colors_R3[condition_percentage + 1]
)
condition_aligned_data[[i]] <- data_run %>%
mutate(Run = paste0("R", i), Run_Color = gradient_color, Run_Linetype = run_linetypes[i])
}
# Combine all aligned data for this condition
aligned_data_combined <- bind_rows(condition_aligned_data)
# Filter data based on max time
aligned_data_combined <- aligned_data_combined %>%
filter(Time_aligned <= max_time_templabel)
# Create the plot
plot <- aligned_data_combined %>%
ggplot(aes(x = Time_aligned, y = Cap_aligned, group = Run, color = Run_Color, linetype = Run_Linetype)) +
geom_line(size = 1.2) +
scale_color_identity(
guide = guide_legend(override.aes = list(color = fixed_colors, linetype = run_linetypes)),
labels = names(fixed_colors)
) +
scale_linetype_identity() +
labs(
x = "Time(s)",
y = "C(fF)",
color = "Runs",
linetype = "Runs"
) +
theme_minimal() +
theme(
axis.text = element_text(size = 20, family = "Times New Roman", face = "bold"),
axis.text.x = element_text(margin = margin(t = 10)),
axis.title.x = element_text(size = 20, family = "Times New Roman", face = "bold", margin = margin(t = 0, r = 40, b = 0, l = 10)),
axis.title.y = element_text(size = 20, family = "Times New Roman", face = "bold", margin = margin(t = 0, r = 40, b = 0, l = 10)),
legend.text = element_text(size = 15, family = "Times New Roman", face = "bold"),
legend.title = element_text(size = 15, family = "Times New Roman", face = "bold"),
legend.position = c(0.7, 0.5),
legend.direction = "vertical",
legend.spacing.y = unit(0.002, 'cm'),
legend.key.size = unit(0.8, 'lines'),
panel.background = element_blank(),
plot.background = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.border = element_rect(color = "black", fill = NA)
) +
scale_x_continuous(
limits = c(0, 850),
breaks = seq(0, 800, by = 400)) +
coord_cartesian(ylim = c(130, 180)) +
scale_y_continuous(breaks = c(130, 150, 170))
plot_list[[condition]] <- plot
}
# Combine the plots into a grid
combined_plot <- cowplot::plot_grid(
plotlist = lapply(plot_list, function(p) {
p
}),
nrow = 11,
ncol = 1,
align = "v"
)
# Save the combined plot
ggsave(
filename = paste0("Combined_Plot_EM", temp_label, ".png"),
plot = combined_plot,
height = 20,
width = 2.5
)
}
# Apply the process_data function to each temperature
process_selected_data(data_25C, "25C", df_25C, df2_25C, max_time_25C)
####################################################
# Function to calculate the maximum cumulative time for the selected run of a temperature
calculate_max_time_for_run <- function(data, df2) {
# Extract the first selected run from df2
selected_run <- str_split(df2$combination[1], ", ")[[1]][1]
# Process the data to calculate cumulative time for the selected run
selected_data <- data %>%
mutate(Condition = str_replace(Condition, "-run", "_run")) %>%
filter(Condition == paste0(df2$condition[1], "_", selected_run)) %>%
group_by(Condition) %>%
mutate(Cumulative_Time = cumsum(Time)) %>%
ungroup()
# Return the maximum cumulative time for the selected run
return(max(selected_data$Cumulative_Time))
}
# Calculate the maximum time for each temperature's selected run
max_time_25C <- calculate_max_time_for_run(data_25C, df2_25C)
# Function to process and visualize data for each temperature using the specific max time
process_selected_data <- function(data, temp_label, df, df2, max_time_templabel) {
data <- data %>%
group_by(Condition) %>%
mutate(Cumulative_Time = cumsum(Time)) %>%
ungroup() %>%
mutate(Time = Cumulative_Time) %>%
select(-Cumulative_Time) %>%
mutate(Condition = str_replace(Condition, "-run", "_run"))
# Ensure 'df' is structured properly
df$condition <- factor(df$condition)
# Preparing to store selected data
selected_data <- data.frame()
# Create ethanol labels for the legend specific to the 25°C condition
Methanol_labels <- c(
"0% Methanol_run5" = "0",
"10% Methanol_run4" = "10",
"20% Methanol_run10" = "20",
"30% Methanol_run5" = "30",
"40% Methanol_run4" = "40",
"50% Methanol_run6" = "50",
"60% Methanol_run5" = "60",
"70% Methanol_run4" = "70",
"80% Methanol_run2" = "80",
"90% Methanol_run3" = "90",
"Pure Methanol_run1" = "100"
)
# Process each condition from df2
unique_conditions <- unique(df2$condition)
for (condition in unique_conditions) {
df2_subset <- df2 %>% filter(condition == !!condition)
selected_runs <- unlist(strsplit(as.character(df2_subset$combination[1]), ", "))
run <- selected_runs[1] # Ensure 'run' exists in our data columns
filtered_data <- data %>%
filter(Condition == paste(condition, run, sep = "_"))
filtered_data <- filtered_data %>%
filter(!is.na(Cap))
# Get the sharpness time for aligning data
sharpness_time <- get_sharpness_time(filtered_data)
data_run <- filtered_data %>%
mutate(Time_aligned = pmax(Time - sharpness_time, 0)) %>%
mutate(Cap_aligned = Cap)
selected_data <- bind_rows(selected_data, data_run %>% select(Time_aligned, Cap_aligned, Condition))
}
# Determine a common starting value, e.g., the median of the first capacitance values of each condition
common_start_value <- median(selected_data$Cap_aligned[selected_data$Time_aligned == 0])
# Function to adjust the capacitance values so that all start from the common start value
adjust_capacitance <- function(data, common_value) {
data %>%
group_by(Condition) %>%
mutate(adjustment = common_value - first(Cap_aligned)) %>%
ungroup() %>%
mutate(Cap_aligned = Cap_aligned + adjustment)
}
# Apply the adjustment function
adjusted_data <- adjust_capacitance(selected_data, common_start_value)
# Truncate data at the max time for the current temperature
adjusted_data <- adjusted_data %>%
filter(Time_aligned <= max_time_templabel)
# Update the Condition column to reflect ethanol percentages
adjusted_data <- adjusted_data %>%
mutate(Condition = Methanol_labels[Condition])
adjusted_data$Condition <- factor(adjusted_data$Condition, levels = c("0", "10", "20", "30", "40", "50", "60", "70", "80", "90", "100"))
# Generate gradient colors from light blue to navy blue for 11 conditions (0% to 100%)
color_gradient <- colorRampPalette(c("#66a61e", "#00441b"))(11)
# Plot the adjusted data with custom color and legend settings
selected_plot <- ggplot(adjusted_data, aes(x = Time_aligned, y = Cap_aligned, color = Condition)) +
geom_line(size = 1) +
labs(
x = "Time(s)",
y = "C(fF)",
color = paste("M_W(%),", temp_label)) + # Adding temp_label to the legend title
theme_minimal() +
theme(
axis.text = element_text(size = 20, family = "Times New Roman", face = "bold"),
axis.text.x = element_text(margin = margin(t = 10)), # Add margin to x-axis text
axis.title.x = element_text(size = 20, family = "Times New Roman", face = "bold", margin = margin(t = 0, r = 10, b = 0, l = 5)),
axis.title.y = element_text(size = 20, family = "Times New Roman", face = "bold", margin = margin(t = 0, r = 10, b = 0, l = 5)),
legend.text = element_text(size = 12, family = "Times New Roman", face = "bold"),
legend.title = element_text(size = 12, family = "Times New Roman", face = "bold"),
legend.position = "none",
legend.direction = "vertical",
legend.spacing.y = unit(0.002, 'cm'),
legend.key.size = unit(0.8, 'lines'),
panel.background = element_blank(),
plot.background = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.border = element_rect(color = "black", fill = NA)
) +
# Apply scale_color_manual with dynamically generated gradient colors
scale_color_manual(
values = color_gradient
)+
scale_x_continuous(
limits = c(0, 850), # Slightly extend the limit
breaks = seq(0, 800, by = 400)) +
coord_cartesian(
ylim = if (temp_label == "60C") c(130, 190) else c(130, 180)
) +
scale_y_continuous(
breaks = if (temp_label == "60C") c(130, 150, 180) else c(130, 150, 170)
)
# Save the plot
ggsave(filename = paste0("Selected_Curves_EM", temp_label, ".png"), plot = selected_plot, height = 2, width = 2.5)
}
# Apply the process_selected_data function to the data of each temperature using its specific max time
process_selected_data(data_25C, "25C", df_25C, df2_25C, max_time_25C)
############
combine_and_align <- function(upper_image_path, lower_image_path, output_path) {
# Load images
upper_image <- image_read(upper_image_path)
lower_image <- image_read(lower_image_path)
# Trim the images to calculate the content width
trimmed_upper_image <- image_trim(upper_image)
trimmed_lower_image <- image_trim(lower_image)
# Get the width of the trimmed images
trimmed_upper_width <- image_info(trimmed_upper_image)$width
trimmed_lower_width <- image_info(trimmed_lower_image)$width
# Calculate the left padding (difference between original width and trimmed width)
left_padding_upper <- image_info(upper_image)$width - trimmed_upper_width
left_padding_lower <- image_info(lower_image)$width - trimmed_lower_width
# Determine the maximum width for alignment
aligned_width <- max(image_info(upper_image)$width, image_info(lower_image)$width)
# Adjust both images with consistent left padding and total width
aligned_upper_image <- image_extent(upper_image, paste0(aligned_width, "x"), gravity = "west")
aligned_lower_image <- image_extent(lower_image, paste0(aligned_width, "x"), gravity = "west")
# Stack images vertically
combined_image <- image_append(c(aligned_upper_image, aligned_lower_image), stack = TRUE)
# Save the final aligned image
image_write(combined_image, path = output_path)
}
# Use the function with your specific files
combine_and_align("Combined_Plot_EM25C.png", "Selected_Curves_EM25C.png", "Combined_Plot_EM25C.png")
#__________________________
#arrange
#______________________________
library(magick)
# Define images with scaling and trimming
images <- list(
image_scale(image_trim(image_read("Combined_Plot_E25C.png")), "2000x3000"),
image_scale(image_trim(image_read("Combined_Plot_E40C.png")), "2000x3000"),
image_scale(image_trim(image_read("Combined_Plot_E50C.png")), "2000x3000"),
image_scale(image_trim(image_read("Combined_Plot_E60C.png")), "2000x3000"),
image_scale(image_trim(image_read("Combined_Plot_M25C.png")), "2000x3000"),
image_scale(image_trim(image_read("Combined_Plot_M40C.png")), "2000x3000"),
image_scale(image_trim(image_read("Combined_Plot_M50C.png")), "2000x3000"),
image_scale(image_trim(image_read("Combined_Plot_M60C.png")), "2000x3000"),
image_scale(image_trim(image_read("Combined_Plot_EM25C.png")), "2000x3000")
)
# Define a transparent spacer (e.g., 50 pixels wide, same height as images)
spacer <- image_blank(width = 15, height = 3000, color = "none")
# Interleave images and spacers
images_with_spacers <- do.call(c, lapply(images, function(img) list(img, spacer)))
# Remove the last spacer (not needed at the end)
images_with_spacers <- images_with_spacers[-length(images_with_spacers)]
# Combine images horizontally
combined_image <- Reduce(function(x, y) image_append(c(x, y), stack = FALSE), images_with_spacers)
# Save the combined image
image_write(combined_image, "Combined_All_Plots.png")
image <- image_read("Combined_All_Plots.png")
# Define the labels
labels <- c("0%", "10%", "20%", "30%", "40%", "50%", "60%", "70%", "80%", "90%", "100%", "All")
# Create a blank label image with the same height as the combined image
label_width <- 95 # Width for the text area
label_height <- image_info(combined_image)$height
label_image <- image_blank(width = label_width, height = label_height, color = "none")
image_height <- image_info(image)$height
image_width <- image_info(image)$width
row_height <- image_height / length(labels)
# Annotate each label in the corresponding position
label_spacing <- label_height / length(labels)
for (i in seq_along(labels)) {
label_y <- (i - 1) * row_height + row_height / 2 # Center align each label
label_image <- image_annotate(
label_image,
text = labels[i],
size = 40, # Adjust font size
gravity = "northwest",
location = paste0("+10+", round(label_y - 55)), # Adjust for fine-tuning vertical centering
color = "black",
font = "Times New Roman", # Font name
weight = 700 # Add this line for bold text (700 is equivalent to bold weight)
)
}
# Combine the label image and the combined image
final_image <- image_append(c(label_image, combined_image), stack = FALSE)
# Save the final image
image_write(final_image, "Combined_All_Plots_Labeled.png")
knitr::include_graphics("Combined_All_Plots_Labeled.png")

1.2 Figure4
rm(list=ls())
#_______________________
# Libraries
#_______________________
library(nortest)
library(writexl)
library(tidyverse)
library(openxlsx)
library(dplyr)
library(Matrix)
library(lme4)
library(ggplot2)
library(car)
library(extrafont)
library(stringr)
library(gridExtra)
library(combinat)
library(tidyr)
library(cowplot)
library(ggpubr)
library(rlang)
library(broom)
library(purrr)
library(ggpmisc) # For adding equation and R² to plot
# Import fonts
# font_import(prompt = FALSE)
loadfonts(device = "win")
#_______________________
# Function Definitions
#_______________________
get_sharpness_time <- function(data) {
# Loop through the data until the cap value increases significantly
for (i in 2:nrow(data)) {
# Check if there's a notable increase in the Cap value after a steady range
if (data$Cap[i] > 160) {
return(data$Time[i])
}
}
# Return NA if no sharpness time is found
return(NA)
}
# Function to extract timeline points and calculate delta times
get_timeline_points <- function(data, condition, run_name) {
# Filter data for the specified condition and run
# Group and calculate cumulative time
data <- data %>%
group_by(Condition) %>%
mutate(Cumulative_Time = cumsum(Time)) %>%
ungroup() %>%
mutate(Time = Cumulative_Time) %>%
select(-Cumulative_Time)
# Replace "-run" with "_run" in the Condition column to ensure uniformity
data <- data %>%
mutate(Condition = str_replace(Condition, "-run", "_run"))
filtered_data <- data %>% filter(Condition == paste(condition, run_name, sep = "_"))
# Get the sharpness time (steepest capacitance transition)
sharpness_time <- get_sharpness_time(filtered_data)
# Filter data based on the sharpness time
data_before_sharpness <- filtered_data %>%
filter(Time < sharpness_time)
# Filter data based on the sharpness time
data_after_sharpness <- filtered_data %>%
filter(Time >= sharpness_time)
# Define timeline points based on the sharpness
max_cap_value <- max(data_after_sharpness$Cap, na.rm = TRUE)
lower_bound <- max_cap_value - 5
upper_bound <- max_cap_value + 5
cap_range_data <- data_after_sharpness %>%
filter(Cap >= lower_bound & Cap <= upper_bound)
first_flat_time <- cap_range_data$Time[1]
last_flat_time <- cap_range_data$Time[nrow(cap_range_data)]
data_after_last_time <- data_after_sharpness %>%
filter(Time > last_flat_time)
delta_t1 <- max(data_before_sharpness$Time) - min(data_before_sharpness$Time)
delta_t2 <- max(cap_range_data$Time) - min(cap_range_data$Time)
delta_t3 <- max(data_after_last_time$Time) - min(data_after_last_time$Time)
TOE <- max(filtered_data$Time)-min(filtered_data$Time)
return(list(
before_flat = data_before_sharpness,
flat = cap_range_data,
after_flat = data_after_last_time,
delta_t1 = delta_t1,
delta_t2 = delta_t2,
delta_t3 = delta_t3,
TOE = TOE
))
}
# Function to calculate Mean and SD of delta_t1, delta_t3, and TOE for each combination of three runs
calculate_mean_sd_combinations <- function(data) {
condition <- unique(data$condition)
# Generate all possible combinations of three runs
run_combinations <- combn(data$run_name, 3, simplify = FALSE)
# Calculate Mean and SD of delta_t1, delta_t3, and TOE for each combination
mean_delta_t1 <- sapply(run_combinations, function(runs) {
selected_delta_t1 <- data %>% filter(run_name %in% runs) %>% pull(delta_t1)
mean(selected_delta_t1)
})
sd_delta_t1 <- sapply(run_combinations, function(runs) {
selected_delta_t1 <- data %>% filter(run_name %in% runs) %>% pull(delta_t1)
sd(selected_delta_t1)
})
mean_delta_t2 <- sapply(run_combinations, function(runs) {
selected_delta_t2 <- data %>% filter(run_name %in% runs) %>% pull(delta_t2)
mean(selected_delta_t2)
})
sd_delta_t2 <- sapply(run_combinations, function(runs) {
selected_delta_t2 <- data %>% filter(run_name %in% runs) %>% pull(delta_t2)
sd(selected_delta_t2)
})
mean_delta_t3 <- sapply(run_combinations, function(runs) {
selected_delta_t3 <- data %>% filter(run_name %in% runs) %>% pull(delta_t3)
mean(selected_delta_t3)
})
sd_delta_t3 <- sapply(run_combinations, function(runs) {
selected_delta_t3 <- data %>% filter(run_name %in% runs) %>% pull(delta_t3)
sd(selected_delta_t3)
})
mean_toe <- sapply(run_combinations, function(runs) {
selected_toes <- data %>% filter(run_name %in% runs) %>% pull(TOE)
mean(selected_toes)
})
sd_toe <- sapply(run_combinations, function(runs) {
selected_toes <- data %>% filter(run_name %in% runs) %>% pull(TOE)
sd(selected_toes)
})
# Create a data frame with the combinations, their Mean and SD for delta_t1, delta_t3, and TOE
mean_sd_table <- data.frame(
condition = condition,
combination = sapply(run_combinations, paste, collapse = ", "),
mean_delta_t1 = mean_delta_t1,
sd_delta_t1 = sd_delta_t1,
mean_delta_t2 = mean_delta_t2,
sd_delta_t2 = sd_delta_t2,
mean_delta_t3 = mean_delta_t3,
sd_delta_t3 = sd_delta_t3,
mean_toe = mean_toe,
sd_toe = sd_toe
)
# Return the Mean and SD table
return(mean_sd_table)
}
# Function to calculate mean Cap for the selected runs
calculate_mean_cap <- function(data, condition, combination) {
# Extract the runs from the combination column
selected_runs <- str_split(combination, ", ") %>% unlist()
# Filter the data for the specific condition and selected runs
filtered_data <- data %>%
filter(str_detect(Condition, paste0(condition, "_(", paste0(selected_runs, collapse = "|"), ")")))
# Calculate the mean Cap
mean_cap <- mean(filtered_data$Cap, na.rm = TRUE)
sd_cap <- sd(filtered_data$Cap, na.rm = TRUE)
return(list(mean_cap = mean_cap, sd_cap = sd_cap))
}
# Function to Process Data for a Given Temperature
process_data <- function(data, temperature) {
# Replace "-run" with "_run" in the Condition column to ensure uniformity
data <- data %>%
mutate(Condition = str_replace(Condition, "-run", "_run"))
unique_conditions <- unique(data$Condition)
split_conditions <- data.frame(
condition = sapply(strsplit(unique_conditions, "_"), `[`, 1),
run_name = sapply(strsplit(unique_conditions, "_"), `[`, 2)
)
results <- split_conditions %>%
mutate(timeline_points = map2(condition, run_name, ~get_timeline_points(data, .x, .y)))
final_results <- results %>%
rowwise() %>%
mutate(
before_flat = list(timeline_points$before_flat),
flat = list(timeline_points$flat),
after_flat = list(timeline_points$after_flat),
delta_t1 = timeline_points$delta_t1,
delta_t2 = timeline_points$delta_t2,
delta_t3 = timeline_points$delta_t3,
TOE = timeline_points$TOE
) %>%
select(-timeline_points) %>%
ungroup()
final_results <- final_results %>%
mutate(
before_flat = sapply(before_flat, function(x) toString(x$Time)),
flat = sapply(flat, function(x) toString(x$Time)),
after_flat = sapply(after_flat, function(x) toString(x$Time))
)
# Save the final results as an Excel file
write_xlsx(final_results, paste0("timeline_results__ethanol_water", temperature, ".xlsx"))
mean_sd_combinations_table <- final_results %>%
group_by(condition) %>%
do(calculate_mean_sd_combinations(.)) %>%
ungroup()
# Find the combination with the smallest combined SD for delta_t1, delta_t2, delta_t3, and TOE
best_combinations_table <- mean_sd_combinations_table %>%
group_by(condition) %>%
mutate(total_sd = sd_delta_t1 + sd_delta_t2 + sd_delta_t3 + sd_toe) %>%
slice_min(total_sd) %>%
ungroup()
# Apply the function to calculate mean Cap for each row of best_combinations_table
best_combinations_table <- best_combinations_table %>%
rowwise() %>%
mutate(
mean_cap = calculate_mean_cap(data, as.character(condition), as.character(combination))$mean_cap,
sd_cap = calculate_mean_cap(data, as.character(condition), as.character(combination))$sd_cap
) %>%
ungroup()
# Define the correct order of conditions
condition_order <- c(
"Pure Water", "10% Ethanol", "20% Ethanol", "30% Ethanol",
"40% Ethanol", "50% Ethanol", "60% Ethanol", "70% Ethanol",
"80% Ethanol", "90% Ethanol", "Pure Ethanol"
)
# Convert the condition column to a factor with the specified levels
best_combinations_table$condition <- factor(best_combinations_table$condition, levels = condition_order)
# Arrange the data frame based on the factor levels
best_combinations_table <- best_combinations_table %>%
arrange(condition)
# Save the tables to Excel
write_xlsx(list(
Mean_SD_Combinations = mean_sd_combinations_table,
Best_Combinations = best_combinations_table
), paste0("mean_sd_combinations_and_selected_runs_ethanol_water", temperature, ".xlsx"))
return(best_combinations_table)
}
# Process for each temperature
data_25C <- read.xlsx("reshape_Water-Ethanol.xlsx", sheet = "25C", colNames = TRUE)
result_25C <- process_data(data_25C, "25C")
##########R2################
condition_order <- c(
"Pure Water", "10% Ethanol", "20% Ethanol", "30% Ethanol",
"40% Ethanol", "50% Ethanol", "60% Ethanol", "70% Ethanol",
"80% Ethanol", "90% Ethanol", "Pure Ethanol"
)
combined_results <- bind_rows(
result_25C %>% mutate(temperature = "25C")
)
plot_cubic_models_25C <- function(data) {
# Filter data for 25C temperature
data_25C <- data %>% filter(temperature == "25C")
# Function to calculate accuracy metrics (R², RMSE, etc.)
calculate_metrics <- function(model, data, y_var) {
predictions <- predict(model, data)
residuals <- data[[y_var]] - predictions
r_squared <- 1 - sum(residuals^2) / sum((data[[y_var]] - mean(data[[y_var]]))^2) # Manually calculate R²
rmse <- sqrt(mean(residuals^2))
mae <- mean(abs(residuals))
range_y <- max(data[[y_var]]) - min(data[[y_var]])
rmse_percentage <- (rmse / range_y) * 100
return(list(r_squared = r_squared, rmse = rmse, mae = mae, rmse_percentage = rmse_percentage))
}
# Fit cubic and LOESS models and calculate metrics for each Delta
lm_delta_t1 <- lm(mean_delta_t1 ~ as.numeric(condition), data = data_25C)
loess_delta_t1 <- loess(mean_delta_t1 ~ as.numeric(condition), data = data_25C, span = 0.75)
metrics_delta_t1_lm <- calculate_metrics(lm_delta_t1, data_25C, "mean_delta_t1")
metrics_delta_t1_loess <- calculate_metrics(loess_delta_t1, data_25C, "mean_delta_t1")
lm_delta_t2 <- lm(mean_delta_t2 ~ as.numeric(condition), data = data_25C)
loess_delta_t2 <- loess(mean_delta_t2 ~ as.numeric(condition), data = data_25C, span = 0.75)
metrics_delta_t2_lm <- calculate_metrics(lm_delta_t2, data_25C, "mean_delta_t2")
metrics_delta_t2_loess <- calculate_metrics(loess_delta_t2, data_25C, "mean_delta_t2")
lm_delta_t3 <- lm(mean_delta_t3 ~ as.numeric(condition), data = data_25C)
loess_delta_t3 <- loess(mean_delta_t3 ~ as.numeric(condition), data = data_25C, span = 0.75)
metrics_delta_t3_lm <- calculate_metrics(lm_delta_t3, data_25C, "mean_delta_t3")
metrics_delta_t3_loess <- calculate_metrics(loess_delta_t3, data_25C, "mean_delta_t3")
lm_toe <- lm(mean_toe ~ as.numeric(condition), data = data_25C)
loess_toe <- loess(mean_toe ~ as.numeric(condition), data = data_25C, span = 0.75)
metrics_toe_lm <- calculate_metrics(lm_toe, data_25C, "mean_toe")
metrics_toe_loess <- calculate_metrics(loess_toe, data_25C, "mean_toe")
# Function to create plots with LOESS and linear model
create_combined_plot <- function(data, y_var, y_label, lm_model, loess_model, metrics_lm, metrics_loess, color) {
min_value <- min(data[[y_var]], na.rm = TRUE)
max_value <- max(data[[y_var]], na.rm = TRUE)
annotation_x <- max(as.numeric(data$condition)) - 1 # Adjust annotation position
values <- if (y_var == "mean_delta_t1") {
list(
max_value = max_value - (max_value - min_value) * 0.1, # Default for mean_delta_t1
legend = c(0.2, 0.8) # Legend position for mean_delta_t1
)
} else {
list(
max_value = max_value - (max_value - min_value) * 0.05, # Move up slightly for other conditions
legend = c(0.2, 0.25) # Legend position for other conditions
)
}
# Access the calculated values
max_value <- values$max_value
legend_position <- values$legend
# Table for R² metrics
table_data <- data.frame(
Model = c("Linear", "LOESS"),
R2 = c(round(metrics_lm$r_squared, 3), round(metrics_loess$r_squared, 3))
)
table_grob <- tableGrob(
table_data,
rows = NULL,
theme = ttheme_minimal(
core = list(bg_params = list(col = "black", lwd = 1),
fg_params = list(fontface = "bold", fontsize = 12)),
colhead = list(bg_params = list(col = "black", lwd = 1),
fg_params = list(fontface = "bold", fontsize = 12))
)
)
ggplot(data, aes(x = as.numeric(condition), y = !!sym(y_var))) +
# Points for "25C"
geom_point(aes(color = "25C"), size = 4, show.legend = TRUE) +
# Linear model (solid line)
# Linear model (solid line)
geom_smooth(method = "lm", formula = y ~ x, se = FALSE,
aes(linetype = "Linear"), color = color, size = 1) +
# LOESS model (dashed line)
geom_line(aes(y = predict(loess_model, data), linetype = "LOESS"),
color = color, size = 1.8) +
# Manual legends and scales
labs(
x = "Ethanol (%)",
y = y_label,
color = "Condition", # Legends for points
linetype = "Model" # Legends for lines
) +
scale_linetype_manual(values = c("Linear" = "solid", "LOESS" = "dashed")) +
scale_color_manual(values = c("25C" = "#B8860B")) +
# Fix the legend guides
guides(
color = guide_legend(override.aes = list(linetype = "blank", shape = 16, size = 4)), # Ensure points are circles
linetype = guide_legend(override.aes = list(color = color, size = c(0.8, 1.5))) # Ensure lines are lines
) +
# Adjust the theme and legend
theme_minimal() +
theme(
axis.text.x = element_text(size = 20, family = "Times New Roman", face = "bold"),
axis.text.y = element_text(size = 20, family = "Times New Roman", face = "bold"),
axis.title = element_text(size = 20, family = "Times New Roman", face = "bold"),
legend.text = element_text(size = 15, family = "Times New Roman", face = "bold"),
legend.title = element_text(size = 15, family = "Times New Roman", face = "bold"),
legend.position = legend_position,
legend.direction = "vertical",
legend.spacing.y = unit(0.002, 'cm'),
#legend.key.size = unit(0.8, 'lines'),
panel.background = element_blank(),
plot.background = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.border = element_rect(color = "black", fill = NA)
)+
annotation_custom(
table_grob,
xmin = max(as.numeric(data$condition)) - 3, # Adjust to position in top-right corner
xmax = max(as.numeric(data$condition))-1.5, # Align to the right edge
ymin = max_value - (max_value - min_value) * 0.1, # Slightly below the top
ymax = max_value # Align to the top edge
)
}
# Combined plots for each metric
p1 <- create_combined_plot(data_25C, "mean_delta_t1", expression(bold(Delta~t[1]~"(s)")), lm_delta_t1, loess_delta_t1, metrics_delta_t1_lm, metrics_delta_t1_loess, "navyblue")
p2 <- create_combined_plot(data_25C, "mean_delta_t2", expression(bold(Delta~t[2]~"(s)")), lm_delta_t2, loess_delta_t2, metrics_delta_t2_lm, metrics_delta_t2_loess, "navyblue")
p3 <- create_combined_plot(data_25C, "mean_delta_t3", expression(bold(Delta~t[3]~"(s)")), lm_delta_t3, loess_delta_t3, metrics_delta_t3_lm, metrics_delta_t3_loess, "navyblue")
p4 <- create_combined_plot(data_25C, "mean_toe", "TOE (s)", lm_toe, loess_toe, metrics_toe_lm, metrics_toe_loess, "navyblue")
plot_grid <- arrangeGrob(p1, p2, p3, p4, ncol = 4)
ggsave("R2E.png", plot = plot_grid, width = 16, height = 6)
}
combined_results$condition <- factor(combined_results$condition, levels = condition_order)
plot_cubic_models_25C(combined_results)
##############RMSE_cap###########
condition_order <- c(
"Pure Water", "10% Ethanol", "20% Ethanol", "30% Ethanol",
"40% Ethanol", "50% Ethanol", "60% Ethanol", "70% Ethanol",
"80% Ethanol", "90% Ethanol", "Pure Ethanol"
)
all_results <- bind_rows(
result_25C %>% mutate(Temperature = "25C")
)
# Now convert the 'Condition' column to a factor with the specified levels
all_results$condition <- factor(all_results$condition, levels = condition_order)
all_results_summarized <- all_results %>%
group_by(condition, Temperature) %>%
summarize(
mean_cap = mean(mean_cap),
sd_cap = mean(sd_cap),
.groups = 'drop' # Ensures we drop the grouping after summarizing
)
# Fit separate models for 25C and 40C using cubic polynomial fit for mean_delta__cap
temp_models_cap <- all_results_summarized %>%
group_by(Temperature) %>%
nest() %>%
mutate(model = map(data, ~ lm(mean_cap ~ as.numeric(condition), data = .x))) # Linear model
# Calculate accuracy metrics for mean_delta__cap
temp_models_cap <- temp_models_cap %>%
mutate(
predictions = map(model, ~ augment(.x)), # Get predictions
r_squared = map_dbl(model, ~ glance(.x)$r.squared), # R-squared
range_y = map_dbl(data, ~ max(.x$mean_cap) - min(.x$mean_cap)), # Range of the observed values
rmse = map_dbl(model, ~ sqrt(mean(residuals(.x)^2))), # Root Mean Squared Error
mae = map_dbl(model, ~ mean(abs(residuals(.x)))) # Mean Absolute Error
) %>%
mutate(
rmse_percentage = (rmse / range_y) * 100, # RMSE as percentage of the range
mae_percentage = (mae / range_y) * 100 # MAE as percentage of the range
)
# Print accuracy metrics for each Temperature for mean_cap
temp_models_cap %>%
select(Temperature, r_squared, rmse, mae, rmse_percentage, mae_percentage)
## # A tibble: 1 × 6
## # Groups: Temperature [1]
## Temperature r_squared rmse mae rmse_percentage mae_percentage
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 25C 0.835 1.94 1.72 12.5 11.1
min_value_cap <- min(all_results_summarized$mean_cap, na.rm = TRUE)
max_value_cap <- max(all_results_summarized$mean_cap, na.rm = TRUE)
max_value_cap = max_value_cap +5
custom_colors <- c("25C" = "navyblue")
# Add LOESS model and compare RMSE with a table
p <- ggplot(all_results_summarized, aes(x = condition, y = mean_cap, group = Temperature)) +
# Points with a fixed color (#B8860B) and a legend
geom_point(aes(color = "25C"), size = 4, show.legend = TRUE) +
# Linear fit with solid navy blue line
geom_smooth(aes(linetype = "Linear"), method = "lm", se = FALSE, size = 1, color = "navyblue") +
# LOESS fit with dashed navy blue line
geom_smooth(aes(linetype = "LOESS"), method = "loess", se = FALSE, size = 1.8, color = "navyblue") +
labs(
x = "Ethanol (%)",
y = expression(bold(Mean~Delta~Cap~(fF))),
color = "Condition", # Title for points legend
linetype = "Model" # Title for lines legend
) +
# Adjust the legend aesthetics
guides(
color = guide_legend(
override.aes = list(
shape = 16, # Circle shape for points
size = 4, # Size of points
color = "#B8860B" # Fixed color for points
)
),
linetype = guide_legend(
override.aes = list(
color = "navyblue", # Navy blue legend lines
size = c(0.8, 1.5) # Thin solid line and thicker dashed line
)
)
) +
# Apply custom scales for linetypes
scale_linetype_manual(values = c("Linear" = "solid", "LOESS" = "dashed")) +
# Apply custom scales for the points legend
scale_color_manual(values = c("25C" = "#B8860B")) +
# Customize x-axis labels
scale_x_discrete(
limits = condition_order,
labels = c("0", "", "20", "", "40", "", "60", "", "80", "", "100") # Show labels every 20
) +
# Apply theme modifications
theme_minimal() +
theme(
axis.text.x = element_text(size = 20, family = "Times New Roman", face = "bold"),
axis.text.y = element_text(size = 20, family = "Times New Roman", face = "bold"),
axis.title = element_text(size = 20, family = "Times New Roman", face = "bold"),
legend.text = element_text(size = 15, family = "Times New Roman", face = "bold"),
legend.title = element_text(size = 15, family = "Times New Roman", face = "bold"),
legend.position = c(0.75, 0.8),
panel.background = element_blank(),
plot.background = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.border = element_rect(color = "black", fill = NA)
)
# Fit LOESS model
loess_model <- loess(mean_cap ~ as.numeric(condition), data = all_results_summarized, span = 0.75)
# Calculate RMSE for LOESS model
loess_predictions <- predict(loess_model, newdata = all_results_summarized)
loess_residuals <- all_results_summarized$mean_cap - loess_predictions
loess_rmse <- sqrt(mean(loess_residuals^2))
loess_range <- max(all_results_summarized$mean_cap) - min(all_results_summarized$mean_cap)
loess_rmse_percentage <- (loess_rmse / loess_range) * 100
# Extract RMSE for linear model
lm_rmse_percentage <- temp_models_cap$rmse_percentage[temp_models_cap$Temperature == "25C"]
# Create RMSE comparison table
rmse_comparison <- data.frame(
Model = c("Linear", "LOESS"),
RMSE = c(paste0(round(lm_rmse_percentage, 3), "%"), paste0(round(loess_rmse_percentage, 3), "%"))
)
#colnames(rmse_comparison) <- c("Model", "RMSE (%)")
# Create table grob
table_grob <- tableGrob(
rmse_comparison,
rows = NULL,
theme = ttheme_minimal(
core = list(bg_params = list(col = "black", lwd = 1),
fg_params = list(fontface = "bold", fontsize = 12)),
colhead = list(bg_params = list(col = "black", lwd = 1),
fg_params = list(fontface = "bold", fontsize = 12))
)
)
# Annotate table in the plot
p <- p +
annotation_custom(
grob = table_grob,
xmin =3, # Position table near the right edge
xmax = 4, # Slightly left of the rightmost edge
ymin = 145, # Slightly below the top
ymax =150 # Align to the top edge
)
# Save the updated plot
ggsave("RMSE_cap_E.png", plot = p, width = 4, height = 6)
#_______________________Metanol_______________________-
# Function to Process Data for a Given Temperature
process_data <- function(data, temperature) {
# Replace "-run" with "_run" in the Condition column to ensure uniformity
data <- data %>%
mutate(Condition = str_replace(Condition, "-run", "_run"))
unique_conditions <- unique(data$Condition)
split_conditions <- data.frame(
condition = sapply(strsplit(unique_conditions, "_"), `[`, 1),
run_name = sapply(strsplit(unique_conditions, "_"), `[`, 2)
)
results <- split_conditions %>%
mutate(timeline_points = map2(condition, run_name, ~get_timeline_points(data, .x, .y)))
final_results <- results %>%
rowwise() %>%
mutate(
before_flat = list(timeline_points$before_flat),
flat = list(timeline_points$flat),
after_flat = list(timeline_points$after_flat),
delta_t1 = timeline_points$delta_t1,
delta_t2 = timeline_points$delta_t2,
delta_t3 = timeline_points$delta_t3,
TOE = timeline_points$TOE
) %>%
select(-timeline_points) %>%
ungroup()
final_results <- final_results %>%
mutate(
before_flat = sapply(before_flat, function(x) toString(x$Time)),
flat = sapply(flat, function(x) toString(x$Time)),
after_flat = sapply(after_flat, function(x) toString(x$Time))
)
# Save the final results as an Excel file
write_xlsx(final_results, paste0("timeline_results__Methanol_water", temperature, ".xlsx"))
mean_sd_combinations_table <- final_results %>%
group_by(condition) %>%
do(calculate_mean_sd_combinations(.)) %>%
ungroup()
# Find the combination with the smallest combined SD for delta_t1, delta_t2, delta_t3, and TOE
best_combinations_table <- mean_sd_combinations_table %>%
group_by(condition) %>%
mutate(total_sd = sd_delta_t1 + sd_delta_t2 + sd_delta_t3 + sd_toe) %>%
slice_min(total_sd) %>%
ungroup()
# Apply the function to calculate mean Cap for each row of best_combinations_table
best_combinations_table <- best_combinations_table %>%
rowwise() %>%
mutate(
mean_cap = calculate_mean_cap(data, as.character(condition), as.character(combination))$mean_cap,
sd_cap = calculate_mean_cap(data, as.character(condition), as.character(combination))$sd_cap
) %>%
ungroup()
# Define the correct order of conditions
condition_order <- c(
"Pure Water", "10% Methanol", "20% Methanol", "30% Methanol",
"40% Methanol", "50% Methanol", "60% Methanol", "70% Methanol",
"80% Methanol", "90% Methanol", "Pure Methanol"
)
# Convert the condition column to a factor with the specified levels
best_combinations_table$condition <- factor(best_combinations_table$condition, levels = condition_order)
# Arrange the data frame based on the factor levels
best_combinations_table <- best_combinations_table %>%
arrange(condition)
# Save the tables to Excel
write_xlsx(list(
Mean_SD_Combinations = mean_sd_combinations_table,
Best_Combinations = best_combinations_table
), paste0("mean_sd_combinations_and_selected_runs_Methanol_water", temperature, ".xlsx"))
return(best_combinations_table)
}
# Process for each temperature
data_25C <- read.xlsx("reshape_Water-Methanol.xlsx", sheet = "25C", colNames = TRUE)
result_25C <- process_data(data_25C, "25C")
##########R2################
condition_order <- c(
"Pure Water", "10% Methanol", "20% Methanol", "30% Methanol",
"40% Methanol", "50% Methanol", "60% Methanol", "70% Methanol",
"80% Methanol", "90% Methanol", "Pure Methanol"
)
combined_results <- bind_rows(
result_25C %>% mutate(temperature = "25C")
)
plot_cubic_models_25C <- function(data) {
# Filter data for 25C temperature
data_25C <- data %>% filter(temperature == "25C")
# Function to calculate accuracy metrics (R², RMSE, etc.)
calculate_metrics <- function(model, data, y_var) {
predictions <- predict(model, data)
residuals <- data[[y_var]] - predictions
r_squared <- 1 - sum(residuals^2) / sum((data[[y_var]] - mean(data[[y_var]]))^2) # Manually calculate R²
rmse <- sqrt(mean(residuals^2))
mae <- mean(abs(residuals))
range_y <- max(data[[y_var]]) - min(data[[y_var]])
rmse_percentage <- (rmse / range_y) * 100
return(list(r_squared = r_squared, rmse = rmse, mae = mae, rmse_percentage = rmse_percentage))
}
# Fit cubic and LOESS models and calculate metrics for each Delta
lm_delta_t1 <- lm(mean_delta_t1 ~ as.numeric(condition), data = data_25C)
loess_delta_t1 <- loess(mean_delta_t1 ~ as.numeric(condition), data = data_25C, span = 0.75)
metrics_delta_t1_lm <- calculate_metrics(lm_delta_t1, data_25C, "mean_delta_t1")
metrics_delta_t1_loess <- calculate_metrics(loess_delta_t1, data_25C, "mean_delta_t1")
lm_delta_t2 <- lm(mean_delta_t2 ~ as.numeric(condition), data = data_25C)
loess_delta_t2 <- loess(mean_delta_t2 ~ as.numeric(condition), data = data_25C, span = 0.75)
metrics_delta_t2_lm <- calculate_metrics(lm_delta_t2, data_25C, "mean_delta_t2")
metrics_delta_t2_loess <- calculate_metrics(loess_delta_t2, data_25C, "mean_delta_t2")
lm_delta_t3 <- lm(mean_delta_t3 ~ as.numeric(condition), data = data_25C)
loess_delta_t3 <- loess(mean_delta_t3 ~ as.numeric(condition), data = data_25C, span = 0.75)
metrics_delta_t3_lm <- calculate_metrics(lm_delta_t3, data_25C, "mean_delta_t3")
metrics_delta_t3_loess <- calculate_metrics(loess_delta_t3, data_25C, "mean_delta_t3")
lm_toe <- lm(mean_toe ~ as.numeric(condition), data = data_25C)
loess_toe <- loess(mean_toe ~ as.numeric(condition), data = data_25C, span = 0.75)
metrics_toe_lm <- calculate_metrics(lm_toe, data_25C, "mean_toe")
metrics_toe_loess <- calculate_metrics(loess_toe, data_25C, "mean_toe")
# Function to create plots with LOESS and linear model
create_combined_plot <- function(data, y_var, y_label, lm_model, loess_model, metrics_lm, metrics_loess, color) {
min_value <- min(data[[y_var]], na.rm = TRUE)
max_value <- max(data[[y_var]], na.rm = TRUE)
annotation_x <- max(as.numeric(data$condition)) - 1 # Adjust annotation position
values <- if (y_var == "mean_delta_t1") {
list(
max_value = max_value - (max_value - min_value) * 0.1, # Default for mean_delta_t1
legend = c(0.2, 0.8) # Legend position for mean_delta_t1
)
} else {
list(
max_value = max_value - (max_value - min_value) * 0.05, # Move up slightly for other conditions
legend = c(0.2, 0.25) # Legend position for other conditions
)
}
# Access the calculated values
max_value <- values$max_value
legend_position <- values$legend
# Table for R² metrics
table_data <- data.frame(
Model = c("Linear", "LOESS"),
R2 = c(round(metrics_lm$r_squared, 3), round(metrics_loess$r_squared, 3))
)
table_grob <- tableGrob(
table_data,
rows = NULL,
theme = ttheme_minimal(
core = list(bg_params = list(col = "black", lwd = 1),
fg_params = list(fontface = "bold", fontsize = 12)),
colhead = list(bg_params = list(col = "black", lwd = 1),
fg_params = list(fontface = "bold", fontsize = 12))
)
)
ggplot(data, aes(x = as.numeric(condition), y = !!sym(y_var))) +
# Points for "25C"
geom_point(aes(color = "25C"), size = 4, show.legend = TRUE) +
# Linear model (solid line)
# Linear model (solid line)
geom_smooth(method = "lm", formula = y ~ x, se = FALSE,
aes(linetype = "Linear"), color = color, size = 1) +
# LOESS model (dashed line)
geom_line(aes(y = predict(loess_model, data), linetype = "LOESS"),
color = color, size = 1.8) +
# Manual legends and scales
labs(
x = "Methanol (%)",
y = y_label,
color = "Condition", # Legends for points
linetype = "Model" # Legends for lines
) +
scale_linetype_manual(values = c("Linear" = "solid", "LOESS" = "dashed")) +
scale_color_manual(values = c("25C" = "#B8860B")) +
# Fix the legend guides
guides(
color = guide_legend(override.aes = list(linetype = "blank", shape = 16, size = 4)), # Ensure points are circles
linetype = guide_legend(override.aes = list(color = color, size = c(0.8, 1.5))) # Ensure lines are lines
) +
# Adjust the theme and legend
theme_minimal() +
theme(
axis.text.x = element_text(size = 20, family = "Times New Roman", face = "bold"),
axis.text.y = element_text(size = 20, family = "Times New Roman", face = "bold"),
axis.title = element_text(size = 20, family = "Times New Roman", face = "bold"),
legend.text = element_text(size = 15, family = "Times New Roman", face = "bold"),
legend.title = element_text(size = 15, family = "Times New Roman", face = "bold"),
legend.position = legend_position,
legend.direction = "vertical",
legend.spacing.y = unit(0.002, 'cm'),
#legend.key.size = unit(0.8, 'lines'),
panel.background = element_blank(),
plot.background = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.border = element_rect(color = "black", fill = NA)
)+
annotation_custom(
table_grob,
xmin = max(as.numeric(data$condition)) - 3, # Adjust to position in top-right corner
xmax = max(as.numeric(data$condition))-1.5, # Align to the right edge
ymin = max_value - (max_value - min_value) * 0.1, # Slightly below the top
ymax = max_value # Align to the top edge
)
}
# Combined plots for each metric
p1 <- create_combined_plot(data_25C, "mean_delta_t1", expression(bold(Delta~t[1]~"(s)")), lm_delta_t1, loess_delta_t1, metrics_delta_t1_lm, metrics_delta_t1_loess, "navyblue")
p2 <- create_combined_plot(data_25C, "mean_delta_t2", expression(bold(Delta~t[2]~"(s)")), lm_delta_t2, loess_delta_t2, metrics_delta_t2_lm, metrics_delta_t2_loess, "navyblue")
p3 <- create_combined_plot(data_25C, "mean_delta_t3", expression(bold(Delta~t[3]~"(s)")), lm_delta_t3, loess_delta_t3, metrics_delta_t3_lm, metrics_delta_t3_loess, "navyblue")
p4 <- create_combined_plot(data_25C, "mean_toe", "TOE (s)", lm_toe, loess_toe, metrics_toe_lm, metrics_toe_loess, "navyblue")
plot_grid <- arrangeGrob(p1, p2, p3, p4, ncol = 4)
ggsave("R2M.png", plot = plot_grid, width = 16, height = 6)
}
combined_results$condition <- factor(combined_results$condition, levels = condition_order)
plot_cubic_models_25C(combined_results)
############rmse_cap#############
###################
all_results <- bind_rows(
result_25C %>% mutate(Temperature = "25C")
)
# Now convert the 'Condition' column to a factor with the specified levels
all_results$condition <- factor(all_results$condition, levels = condition_order)
all_results_summarized <- all_results %>%
group_by(condition, Temperature) %>%
summarize(
mean_cap = mean(mean_cap),
sd_cap = mean(sd_cap),
.groups = 'drop' # Ensures we drop the grouping after summarizing
)
# Fit separate models for 25C and 40C using cubic polynomial fit for mean_delta__cap
temp_models_cap <- all_results_summarized %>%
group_by(Temperature) %>%
nest() %>%
mutate(model = map(data, ~ lm(mean_cap ~ as.numeric(condition), data = .x))) # Linear model
# Calculate accuracy metrics for mean_delta__cap
temp_models_cap <- temp_models_cap %>%
mutate(
predictions = map(model, ~ augment(.x)), # Get predictions
r_squared = map_dbl(model, ~ glance(.x)$r.squared), # R-squared
range_y = map_dbl(data, ~ max(.x$mean_cap) - min(.x$mean_cap)), # Range of the observed values
rmse = map_dbl(model, ~ sqrt(mean(residuals(.x)^2))), # Root Mean Squared Error
mae = map_dbl(model, ~ mean(abs(residuals(.x)))) # Mean Absolute Error
) %>%
mutate(
rmse_percentage = (rmse / range_y) * 100, # RMSE as percentage of the range
mae_percentage = (mae / range_y) * 100 # MAE as percentage of the range
)
# Print accuracy metrics for each Temperature for mean_cap
temp_models_cap %>%
select(Temperature, r_squared, rmse, mae, rmse_percentage, mae_percentage)
## # A tibble: 1 × 6
## # Groups: Temperature [1]
## Temperature r_squared rmse mae rmse_percentage mae_percentage
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 25C 0.374 4.58 3.65 21.3 16.9
min_value_cap <- min(all_results_summarized$mean_cap, na.rm = TRUE)
max_value_cap <- max(all_results_summarized$mean_cap, na.rm = TRUE)
max_value_cap = max_value_cap +5
custom_colors <- c("25C" = "navyblue")
p <- ggplot(all_results_summarized, aes(x = condition, y = mean_cap, group = Temperature)) +
# Points with a fixed color (#B8860B) and a legend
geom_point(aes(color = "25C"), size = 4, show.legend = TRUE) +
# Linear fit with solid navy blue line
geom_smooth(aes(linetype = "Linear"), method = "lm", se = FALSE, size = 1, color = "navyblue") +
# LOESS fit with dashed navy blue line
geom_smooth(aes(linetype = "LOESS"), method = "loess", se = FALSE, size = 1.8, color = "navyblue") +
labs(
x = "Methanol (%)",
y = expression(bold(Mean~Delta~Cap~(fF))),
color = "Condition", # Title for points legend
linetype = "Model" # Title for lines legend
) +
# Adjust the legend aesthetics
guides(
color = guide_legend(
override.aes = list(
shape = 16, # Circle shape for points
size = 4, # Size of points
color = "#B8860B" # Fixed color for points
)
),
linetype = guide_legend(
override.aes = list(
color = "navyblue", # Navy blue legend lines
size = c(0.8, 1.5) # Thin solid line and thicker dashed line
)
)
) +
# Apply custom scales for linetypes
scale_linetype_manual(values = c("Linear" = "solid", "LOESS" = "dashed")) +
# Apply custom scales for the points legend
scale_color_manual(values = c("25C" = "#B8860B")) +
# Customize x-axis labels
scale_x_discrete(
limits = condition_order,
labels = c("0", "", "20", "", "40", "", "60", "", "80", "", "100") # Show labels every 20
) +
# Apply theme modifications
theme_minimal() +
theme(
axis.text.x = element_text(size = 20, family = "Times New Roman", face = "bold"),
axis.text.y = element_text(size = 20, family = "Times New Roman", face = "bold"),
axis.title = element_text(size = 20, family = "Times New Roman", face = "bold"),
legend.text = element_text(size = 15, family = "Times New Roman", face = "bold"),
legend.title = element_text(size = 15, family = "Times New Roman", face = "bold"),
legend.position = c(0.75, 0.8),
panel.background = element_blank(),
plot.background = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.border = element_rect(color = "black", fill = NA)
)
# Fit LOESS model
loess_model <- loess(mean_cap ~ as.numeric(condition), data = all_results_summarized, span = 0.75)
# Calculate RMSE for LOESS model
loess_predictions <- predict(loess_model, newdata = all_results_summarized)
loess_residuals <- all_results_summarized$mean_cap - loess_predictions
loess_rmse <- sqrt(mean(loess_residuals^2))
loess_range <- max(all_results_summarized$mean_cap) - min(all_results_summarized$mean_cap)
loess_rmse_percentage <- (loess_rmse / loess_range) * 100
# Extract RMSE for linear model
lm_rmse_percentage <- temp_models_cap$rmse_percentage[temp_models_cap$Temperature == "25C"]
# Create RMSE comparison table
rmse_comparison <- data.frame(
Model = c("Linear", "LOESS"),
RMSE = c(paste0(round(lm_rmse_percentage, 3), "%"), paste0(round(loess_rmse_percentage, 3), "%"))
)
#colnames(rmse_comparison) <- c("Model", "RMSE (%)")
# Create table grob
table_grob <- tableGrob(
rmse_comparison,
rows = NULL,
theme = ttheme_minimal(
core = list(bg_params = list(col = "black", lwd = 1),
fg_params = list(fontface = "bold", fontsize = 12)),
colhead = list(bg_params = list(col = "black", lwd = 1),
fg_params = list(fontface = "bold", fontsize = 12))
)
)
# Annotate table in the plot
p <- p +
annotation_custom(
grob = table_grob,
xmin =3, # Position table near the right edge
xmax = 4, # Slightly left of the rightmost edge
ymin = 148, # Slightly below the top
ymax =150 # Align to the top edge
)
ggsave("RMSE_cap_M.png", plot = p, width =4, height = 6)
#___________________M_E____________________
# Function to Process Data for a Given Temperature
process_data <- function(data, temperature) {
# Replace "-run" with "_run" in the Condition column to ensure uniformity
data <- data %>%
mutate(Condition = str_replace(Condition, "-run", "_run"))
unique_conditions <- unique(data$Condition)
split_conditions <- data.frame(
condition = sapply(strsplit(unique_conditions, "_"), `[`, 1),
run_name = sapply(strsplit(unique_conditions, "_"), `[`, 2)
)
results <- split_conditions %>%
mutate(timeline_points = map2(condition, run_name, ~get_timeline_points(data, .x, .y)))
final_results <- results %>%
rowwise() %>%
mutate(
before_flat = list(timeline_points$before_flat),
flat = list(timeline_points$flat),
after_flat = list(timeline_points$after_flat),
delta_t1 = timeline_points$delta_t1,
delta_t2 = timeline_points$delta_t2,
delta_t3 = timeline_points$delta_t3,
TOE = timeline_points$TOE
) %>%
select(-timeline_points) %>%
ungroup()
final_results <- final_results %>%
mutate(
before_flat = sapply(before_flat, function(x) toString(x$Time)),
flat = sapply(flat, function(x) toString(x$Time)),
after_flat = sapply(after_flat, function(x) toString(x$Time))
)
# Save the final results as an Excel file
write_xlsx(final_results, paste0("timeline_results__Methanol_Ethanol", temperature, ".xlsx"))
mean_sd_combinations_table <- final_results %>%
group_by(condition) %>%
do(calculate_mean_sd_combinations(.)) %>%
ungroup()
# Find the combination with the smallest combined SD for delta_t1, delta_t2, delta_t3, and TOE
best_combinations_table <- mean_sd_combinations_table %>%
group_by(condition) %>%
mutate(total_sd = sd_delta_t1 + sd_delta_t2 + sd_delta_t3 + sd_toe) %>%
slice_min(total_sd) %>%
ungroup()
# Apply the function to calculate mean Cap for each row of best_combinations_table
best_combinations_table <- best_combinations_table %>%
rowwise() %>%
mutate(
mean_cap = calculate_mean_cap(data, as.character(condition), as.character(combination))$mean_cap,
sd_cap = calculate_mean_cap(data, as.character(condition), as.character(combination))$sd_cap
) %>%
ungroup()
# Define the correct order of conditions
condition_order <- c(
"0% Methanol", "10% Methanol", "20% Methanol", "30% Methanol",
"40% Methanol", "50% Methanol", "60% Methanol", "70% Methanol",
"80% Methanol", "90% Methanol", "Pure Methanol"
)
# Convert the condition column to a factor with the specified levels
best_combinations_table$condition <- factor(best_combinations_table$condition, levels = condition_order)
# Arrange the data frame based on the factor levels
best_combinations_table <- best_combinations_table %>%
arrange(condition)
# Save the tables to Excel
write_xlsx(list(
Mean_SD_Combinations = mean_sd_combinations_table,
Best_Combinations = best_combinations_table
), paste0("mean_sd_combinations_and_selected_runs_Methanol_Ethanol", temperature, ".xlsx"))
return(best_combinations_table)
}
# Process for each temperature
data_25C <- read.xlsx("reshape_Methanol-Ethanol.xlsx", sheet = "25C", colNames = TRUE)
result_25C <- process_data(data_25C, "25C")
##############R2##################
combined_results <- bind_rows(
result_25C %>% mutate(temperature = "25C")
)
all_results_summarized <- combined_results %>%
group_by(condition, temperature) %>%
summarize(
mean_cap = mean(mean_cap),
sd_cap = mean(sd_cap),
.groups = 'drop' # Ensures we drop the grouping after summarizing
)
min_value_cap <- min(all_results_summarized$mean_cap, na.rm = TRUE)
max_value_cap <- max(all_results_summarized$mean_cap, na.rm = TRUE)
plot_cubic_models_25C <- function(data) {
# Filter data for 25C temperature
data_25C <- data %>% filter(temperature == "25C")
# Function to calculate accuracy metrics (R², RMSE, etc.)
calculate_metrics <- function(model, data, y_var) {
predictions <- predict(model, data)
residuals <- data[[y_var]] - predictions
r_squared <- 1 - sum(residuals^2) / sum((data[[y_var]] - mean(data[[y_var]]))^2) # Manually calculate R²
rmse <- sqrt(mean(residuals^2))
mae <- mean(abs(residuals))
range_y <- max(data[[y_var]]) - min(data[[y_var]])
rmse_percentage <- (rmse / range_y) * 100
return(list(r_squared = r_squared, rmse = rmse, mae = mae, rmse_percentage = rmse_percentage))
}
# Fit cubic and LOESS models and calculate metrics for each Delta
lm_delta_t1 <- lm(mean_delta_t1 ~ as.numeric(condition), data = data_25C)
loess_delta_t1 <- loess(mean_delta_t1 ~ as.numeric(condition), data = data_25C, span = 0.75)
metrics_delta_t1_lm <- calculate_metrics(lm_delta_t1, data_25C, "mean_delta_t1")
metrics_delta_t1_loess <- calculate_metrics(loess_delta_t1, data_25C, "mean_delta_t1")
lm_delta_t2 <- lm(mean_delta_t2 ~ as.numeric(condition), data = data_25C)
loess_delta_t2 <- loess(mean_delta_t2 ~ as.numeric(condition), data = data_25C, span = 0.75)
metrics_delta_t2_lm <- calculate_metrics(lm_delta_t2, data_25C, "mean_delta_t2")
metrics_delta_t2_loess <- calculate_metrics(loess_delta_t2, data_25C, "mean_delta_t2")
lm_delta_t3 <- lm(mean_delta_t3 ~ as.numeric(condition), data = data_25C)
loess_delta_t3 <- loess(mean_delta_t3 ~ as.numeric(condition), data = data_25C, span = 0.75)
metrics_delta_t3_lm <- calculate_metrics(lm_delta_t3, data_25C, "mean_delta_t3")
metrics_delta_t3_loess <- calculate_metrics(loess_delta_t3, data_25C, "mean_delta_t3")
lm_toe <- lm(mean_toe ~ as.numeric(condition), data = data_25C)
loess_toe <- loess(mean_toe ~ as.numeric(condition), data = data_25C, span = 0.75)
metrics_toe_lm <- calculate_metrics(lm_toe, data_25C, "mean_toe")
metrics_toe_loess <- calculate_metrics(loess_toe, data_25C, "mean_toe")
# Function to create plots with LOESS and linear model
create_combined_plot <- function(data, y_var, y_label, lm_model, loess_model, metrics_lm, metrics_loess, color) {
min_value <- min(data[[y_var]], na.rm = TRUE)
max_value <- max(data[[y_var]], na.rm = TRUE)
annotation_x <- max(as.numeric(data$condition)) - 1 # Adjust annotation position
values <- if (y_var == "mean_delta_t1") {
list(
max_value = max_value - (max_value - min_value) * 0.1, # Default for mean_delta_t1
legend = c(0.2, 0.8) # Legend position for mean_delta_t1
)
} else {
list(
max_value = max_value - (max_value - min_value) * 0.05, # Move up slightly for other conditions
legend = c(0.8, 0.6) # Legend position for other conditions
)
}
# Access the calculated values
max_value <- values$max_value
legend_position <- values$legend
# Table for R² metrics
table_data <- data.frame(
Model = c("Linear", "LOESS"),
R2 = c(round(metrics_lm$r_squared, 3), round(metrics_loess$r_squared, 3))
)
table_grob <- tableGrob(
table_data,
rows = NULL,
theme = ttheme_minimal(
core = list(bg_params = list(col = "black", lwd = 1),
fg_params = list(fontface = "bold", fontsize = 12)),
colhead = list(bg_params = list(col = "black", lwd = 1),
fg_params = list(fontface = "bold", fontsize = 12))
)
)
ggplot(data, aes(x = as.numeric(condition), y = !!sym(y_var))) +
# Points for "25C"
geom_point(aes(color = "25C"), size = 4, show.legend = TRUE) +
# Linear model (solid line)
# Linear model (solid line)
geom_smooth(method = "lm", formula = y ~ x, se = FALSE,
aes(linetype = "Linear"), color = color, size = 1) +
# LOESS model (dashed line)
geom_line(aes(y = predict(loess_model, data), linetype = "LOESS"),
color = color, size = 1.8) +
# Manual legends and scales
labs(
x = "Methanol (%)",
y = y_label,
color = "Condition", # Legends for points
linetype = "Model" # Legends for lines
) +
scale_linetype_manual(values = c("Linear" = "solid", "LOESS" = "dashed")) +
scale_color_manual(values = c("25C" = "#B8860B")) +
# Fix the legend guides
guides(
color = guide_legend(override.aes = list(linetype = "blank", shape = 16, size = 4)), # Ensure points are circles
linetype = guide_legend(override.aes = list(color = color, size = c(0.8, 1.5))) # Ensure lines are lines
) +
# Adjust the theme and legend
theme_minimal() +
theme(
axis.text.x = element_text(size = 20, family = "Times New Roman", face = "bold"),
axis.text.y = element_text(size = 20, family = "Times New Roman", face = "bold"),
axis.title = element_text(size = 20, family = "Times New Roman", face = "bold"),
legend.text = element_text(size = 15, family = "Times New Roman", face = "bold"),
legend.title = element_text(size = 15, family = "Times New Roman", face = "bold"),
legend.position = legend_position,
legend.direction = "vertical",
legend.spacing.y = unit(0.002, 'cm'),
#legend.key.size = unit(0.8, 'lines'),
panel.background = element_blank(),
plot.background = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.border = element_rect(color = "black", fill = NA)
)+
annotation_custom(
table_grob,
xmin = max(as.numeric(data$condition)) - 3, # Adjust to position in top-right corner
xmax = max(as.numeric(data$condition))-1.5, # Align to the right edge
ymin = max_value - (max_value - min_value) * 0.1, # Slightly below the top
ymax = max_value # Align to the top edge
)
}
# Combined plots for each metric
p1 <- create_combined_plot(data_25C, "mean_delta_t1", expression(bold(Delta~t[1]~"(s)")), lm_delta_t1, loess_delta_t1, metrics_delta_t1_lm, metrics_delta_t1_loess, "navyblue")
p2 <- create_combined_plot(data_25C, "mean_delta_t2", expression(bold(Delta~t[2]~"(s)")), lm_delta_t2, loess_delta_t2, metrics_delta_t2_lm, metrics_delta_t2_loess, "navyblue")
p3 <- create_combined_plot(data_25C, "mean_delta_t3", expression(bold(Delta~t[3]~"(s)")), lm_delta_t3, loess_delta_t3, metrics_delta_t3_lm, metrics_delta_t3_loess, "navyblue")
p4 <- create_combined_plot(data_25C, "mean_toe", "TOE (s)", lm_toe, loess_toe, metrics_toe_lm, metrics_toe_loess, "navyblue")
plot_grid <- arrangeGrob(p1, p2, p3, p4, ncol = 4)
ggsave("R2EM.png", plot = plot_grid, width = 16, height = 6)
}
condition_order <- c(
"0% Methanol", "10% Methanol", "20% Methanol", "30% Methanol",
"40% Methanol", "50% Methanol", "60% Methanol", "70% Methanol",
"80% Methanol", "90% Methanol", "Pure Methanol"
)
combined_results$condition <- factor(combined_results$condition, levels = condition_order)
plot_cubic_models_25C(combined_results)
#######################RMSE_CAP###################
condition_order <- c(
"0% Methanol", "10% Methanol", "20% Methanol", "30% Methanol",
"40% Methanol", "50% Methanol", "60% Methanol", "70% Methanol",
"80% Methanol", "90% Methanol", "Pure Methanol"
)
all_results <- bind_rows(
result_25C %>% mutate(Temperature = "25C")
)
# Now convert the 'Condition' column to a factor with the specified levels
all_results$condition <- factor(all_results$condition, levels = condition_order)
all_results_summarized <- all_results %>%
group_by(condition, Temperature) %>%
summarize(
mean_cap = mean(mean_cap),
sd_cap = mean(sd_cap),
.groups = 'drop' # Ensures we drop the grouping after summarizing
)
# Fit separate models for 25C and 40C using cubic polynomial fit for mean_delta__cap
temp_models_cap <- all_results_summarized %>%
group_by(Temperature) %>%
nest() %>%
mutate(model = map(data, ~ lm(mean_cap ~ as.numeric(condition), data = .x))) # Linear model
# Calculate accuracy metrics for mean_delta__cap
temp_models_cap <- temp_models_cap %>%
mutate(
predictions = map(model, ~ augment(.x)), # Get predictions
r_squared = map_dbl(model, ~ glance(.x)$r.squared), # R-squared
range_y = map_dbl(data, ~ max(.x$mean_cap) - min(.x$mean_cap)), # Range of the observed values
rmse = map_dbl(model, ~ sqrt(mean(residuals(.x)^2))), # Root Mean Squared Error
mae = map_dbl(model, ~ mean(abs(residuals(.x)))) # Mean Absolute Error
) %>%
mutate(
rmse_percentage = (rmse / range_y) * 100, # RMSE as percentage of the range
mae_percentage = (mae / range_y) * 100 # MAE as percentage of the range
)
# Print accuracy metrics for each Temperature for mean_cap
temp_models_cap %>%
select(Temperature, r_squared, rmse, mae, rmse_percentage, mae_percentage)
## # A tibble: 1 × 6
## # Groups: Temperature [1]
## Temperature r_squared rmse mae rmse_percentage mae_percentage
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 25C 0.429 1.76 1.38 18.6 14.6
min_value_cap <- min(all_results_summarized$mean_cap, na.rm = TRUE)
max_value_cap <- max(all_results_summarized$mean_cap, na.rm = TRUE)
max_value_cap = max_value_cap +5
custom_colors <- c("25C" = "navyblue")
# Add LOESS model and compare RMSE with a table
p <- ggplot(all_results_summarized, aes(x = condition, y = mean_cap, group = Temperature)) +
# Points with a fixed color (#B8860B) and a legend
geom_point(aes(color = "25C"), size = 4, show.legend = TRUE) +
# Linear fit with solid navy blue line
geom_smooth(aes(linetype = "Linear"), method = "lm", se = FALSE, size = 1, color = "navyblue") +
# LOESS fit with dashed navy blue line
geom_smooth(aes(linetype = "LOESS"), method = "loess", se = FALSE, size = 1.8, color = "navyblue") +
labs(
x = "Methanol (%)",
y = expression(bold(Mean~Delta~Cap~(fF))),
color = "Condition", # Title for points legend
linetype = "Model" # Title for lines legend
) +
# Adjust the legend aesthetics
guides(
color = guide_legend(
override.aes = list(
shape = 16, # Circle shape for points
size = 4, # Size of points
color = "#B8860B" # Fixed color for points
)
),
linetype = guide_legend(
override.aes = list(
color = "navyblue", # Navy blue legend lines
size = c(0.8, 1.5) # Thin solid line and thicker dashed line
)
)
) +
# Apply custom scales for linetypes
scale_linetype_manual(values = c("Linear" = "solid", "LOESS" = "dashed")) +
# Apply custom scales for the points legend
scale_color_manual(values = c("25C" = "#B8860B")) +
# Customize x-axis labels
scale_x_discrete(
limits = condition_order,
labels = c("0", "", "20", "", "40", "", "60", "", "80", "", "100") # Show labels every 20
) +
# Apply theme modifications
theme_minimal() +
theme(
axis.text.x = element_text(size = 20, family = "Times New Roman", face = "bold"),
axis.text.y = element_text(size = 20, family = "Times New Roman", face = "bold"),
axis.title = element_text(size = 20, family = "Times New Roman", face = "bold"),
legend.text = element_text(size = 15, family = "Times New Roman", face = "bold"),
legend.title = element_text(size = 15, family = "Times New Roman", face = "bold"),
legend.position = c(0.75, 0.8),
panel.background = element_blank(),
plot.background = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.border = element_rect(color = "black", fill = NA)
)
# Fit LOESS model
loess_model <- loess(mean_cap ~ as.numeric(condition), data = all_results_summarized, span = 0.75)
# Calculate RMSE for LOESS model
loess_predictions <- predict(loess_model, newdata = all_results_summarized)
loess_residuals <- all_results_summarized$mean_cap - loess_predictions
loess_rmse <- sqrt(mean(loess_residuals^2))
loess_range <- max(all_results_summarized$mean_cap) - min(all_results_summarized$mean_cap)
loess_rmse_percentage <- (loess_rmse / loess_range) * 100
# Extract RMSE for linear model
lm_rmse_percentage <- temp_models_cap$rmse_percentage[temp_models_cap$Temperature == "25C"]
# Create RMSE comparison table
rmse_comparison <- data.frame(
Model = c("Linear", "LOESS"),
RMSE = c(paste0(round(lm_rmse_percentage, 3), "%"), paste0(round(loess_rmse_percentage, 3), "%"))
)
#colnames(rmse_comparison) <- c("Model", "RMSE (%)")
table_grob <- tableGrob(
rmse_comparison,
rows = NULL,
theme = ttheme_minimal(
core = list(bg_params = list(col = "black", lwd = 1),
fg_params = list(fontface = "bold", fontsize = 12)),
colhead = list(bg_params = list(col = "black", lwd = 1),
fg_params = list(fontface = "bold", fontsize = 12))
)
)
# Annotate table in the plot
p <- p +
annotation_custom(
grob = table_grob,
xmin =3.5, # Position table near the right edge
xmax = 4, # Slightly left of the rightmost edge
ymin = 143, # Slightly below the top
ymax =145 # Align to the top edge
)
# Save the plot
ggsave("RMSE_cap_EM.png", plot = p, width = 4, height = 6)
#_______________________-arrange
#___________________________________
library(magick)
# Load and resize the images
R2E <- image_read("R2E.png") %>% image_background("white")
RMSE_cap_E <- image_read("RMSE_cap_E.png") %>% image_background("white")
R2M <- image_read("R2M.png") %>% image_background("white")
RMSE_cap_M <- image_read("RMSE_cap_M.png") %>% image_background("white")
R2EM <- image_read("R2EM.png") %>% image_background("white")
RMSE_cap_EM <- image_read("RMSE_cap_EM.png") %>% image_background("white")
# Create horizontal spacer
horizontal_spacer <- image_blank(width = 100, height = image_info(R2E)$height, color = "white") # 50px wide blank space
# Arrange the first row with horizontal spacers: R2E, spacer, RMSE_cap_E, spacer, RMSE_t2_E
row1 <- image_append(c(R2E, horizontal_spacer, RMSE_cap_E))
# Arrange the second row with horizontal spacers: R2M, spacer, RMSE_cap_M, spacer, RMSE_t2_M
row2 <- image_append(c(R2M, horizontal_spacer, RMSE_cap_M))
# Arrange the third row with horizontal spacers: R2EM, spacer, RMSE_cap_EM, spacer, RMSE_t2_EM
row3 <- image_append(c(R2EM, horizontal_spacer, RMSE_cap_EM))
# Create vertical spacer for alignment
vertical_spacer <- image_blank(width = image_info(row1)$width, height = 200, color = "white") # 200px height blank space
# Stack the rows vertically with spacers between them
final_image <- image_append(c(row1, vertical_spacer, row2, vertical_spacer, row3), stack = TRUE)
# Save the final combined image
image_write(final_image, "combined_plot2.png")
knitr::include_graphics("combined_plot2.png")

1.3 Figure5
rm(list=ls())
#_______________________
# Libraries
#_______________________
library(nortest)
library(writexl)
library(tidyverse)
library(openxlsx)
library(dplyr)
library(Matrix)
library(lme4)
library(ggplot2)
library(car)
library(extrafont)
library(stringr)
library(gridExtra)
library(combinat)
library(tidyr)
library(cowplot)
library(ggpubr)
library(rlang)
library(broom)
library(purrr)
library(ggpmisc) # For adding equation and R² to plot
# Import fonts
# font_import(prompt = FALSE)
loadfonts(device = "win")
#_______________________
# Function Definitions
#_______________________
get_sharpness_time <- function(data) {
# Loop through the data until the cap value increases significantly
for (i in 2:nrow(data)) {
# Check if there's a notable increase in the Cap value after a steady range
if (data$Cap[i] > 160) {
return(data$Time[i])
}
}
# Return NA if no sharpness time is found
return(NA)
}
# Function to extract timeline points and calculate delta times
get_timeline_points <- function(data, condition, run_name) {
# Filter data for the specified condition and run
# Group and calculate cumulative time
data <- data %>%
group_by(Condition) %>%
mutate(Cumulative_Time = cumsum(Time)) %>%
ungroup() %>%
mutate(Time = Cumulative_Time) %>%
select(-Cumulative_Time)
# Replace "-run" with "_run" in the Condition column to ensure uniformity
data <- data %>%
mutate(Condition = str_replace(Condition, "-run", "_run"))
filtered_data <- data %>% filter(Condition == paste(condition, run_name, sep = "_"))
# Get the sharpness time (steepest capacitance transition)
sharpness_time <- get_sharpness_time(filtered_data)
# Filter data based on the sharpness time
data_before_sharpness <- filtered_data %>%
filter(Time < sharpness_time)
# Filter data based on the sharpness time
data_after_sharpness <- filtered_data %>%
filter(Time >= sharpness_time)
# Define timeline points based on the sharpness
max_cap_value <- max(data_after_sharpness$Cap, na.rm = TRUE)
lower_bound <- max_cap_value - 5
upper_bound <- max_cap_value + 5
cap_range_data <- data_after_sharpness %>%
filter(Cap >= lower_bound & Cap <= upper_bound)
first_flat_time <- cap_range_data$Time[1]
last_flat_time <- cap_range_data$Time[nrow(cap_range_data)]
data_after_last_time <- data_after_sharpness %>%
filter(Time > last_flat_time)
delta_t1 <- max(data_before_sharpness$Time) - min(data_before_sharpness$Time)
delta_t2 <- max(cap_range_data$Time) - min(cap_range_data$Time)
delta_t3 <- max(data_after_last_time$Time) - min(data_after_last_time$Time)
TOE <- max(filtered_data$Time)-min(filtered_data$Time)
return(list(
before_flat = data_before_sharpness,
flat = cap_range_data,
after_flat = data_after_last_time,
delta_t1 = delta_t1,
delta_t2 = delta_t2,
delta_t3 = delta_t3,
TOE = TOE
))
}
# Function to calculate Mean and SD of delta_t1, delta_t3, and TOE for each combination of three runs
calculate_mean_sd_combinations <- function(data) {
condition <- unique(data$condition)
# Generate all possible combinations of three runs
run_combinations <- combn(data$run_name, 3, simplify = FALSE)
# Calculate Mean and SD of delta_t1, delta_t3, and TOE for each combination
mean_delta_t1 <- sapply(run_combinations, function(runs) {
selected_delta_t1 <- data %>% filter(run_name %in% runs) %>% pull(delta_t1)
mean(selected_delta_t1)
})
sd_delta_t1 <- sapply(run_combinations, function(runs) {
selected_delta_t1 <- data %>% filter(run_name %in% runs) %>% pull(delta_t1)
sd(selected_delta_t1)
})
mean_delta_t2 <- sapply(run_combinations, function(runs) {
selected_delta_t2 <- data %>% filter(run_name %in% runs) %>% pull(delta_t2)
mean(selected_delta_t2)
})
sd_delta_t2 <- sapply(run_combinations, function(runs) {
selected_delta_t2 <- data %>% filter(run_name %in% runs) %>% pull(delta_t2)
sd(selected_delta_t2)
})
mean_delta_t3 <- sapply(run_combinations, function(runs) {
selected_delta_t3 <- data %>% filter(run_name %in% runs) %>% pull(delta_t3)
mean(selected_delta_t3)
})
sd_delta_t3 <- sapply(run_combinations, function(runs) {
selected_delta_t3 <- data %>% filter(run_name %in% runs) %>% pull(delta_t3)
sd(selected_delta_t3)
})
mean_toe <- sapply(run_combinations, function(runs) {
selected_toes <- data %>% filter(run_name %in% runs) %>% pull(TOE)
mean(selected_toes)
})
sd_toe <- sapply(run_combinations, function(runs) {
selected_toes <- data %>% filter(run_name %in% runs) %>% pull(TOE)
sd(selected_toes)
})
# Create a data frame with the combinations, their Mean and SD for delta_t1, delta_t3, and TOE
mean_sd_table <- data.frame(
condition = condition,
combination = sapply(run_combinations, paste, collapse = ", "),
mean_delta_t1 = mean_delta_t1,
sd_delta_t1 = sd_delta_t1,
mean_delta_t2 = mean_delta_t2,
sd_delta_t2 = sd_delta_t2,
mean_delta_t3 = mean_delta_t3,
sd_delta_t3 = sd_delta_t3,
mean_toe = mean_toe,
sd_toe = sd_toe
)
# Return the Mean and SD table
return(mean_sd_table)
}
# Function to calculate mean Cap for the selected runs
calculate_mean_cap <- function(data, condition, combination) {
# Extract the runs from the combination column
selected_runs <- str_split(combination, ", ") %>% unlist()
# Filter the data for the specific condition and selected runs
filtered_data <- data %>%
filter(str_detect(Condition, paste0(condition, "_(", paste0(selected_runs, collapse = "|"), ")")))
# Calculate the mean Cap
mean_cap <- mean(filtered_data$Cap, na.rm = TRUE)
sd_cap <- sd(filtered_data$Cap, na.rm = TRUE)
return(list(mean_cap = mean_cap, sd_cap = sd_cap))
}
# Function to Process Data for a Given Temperature
process_data <- function(data, temperature) {
# Replace "-run" with "_run" in the Condition column to ensure uniformity
data <- data %>%
mutate(Condition = str_replace(Condition, "-run", "_run"))
unique_conditions <- unique(data$Condition)
split_conditions <- data.frame(
condition = sapply(strsplit(unique_conditions, "_"), `[`, 1),
run_name = sapply(strsplit(unique_conditions, "_"), `[`, 2)
)
results <- split_conditions %>%
mutate(timeline_points = map2(condition, run_name, ~get_timeline_points(data, .x, .y)))
final_results <- results %>%
rowwise() %>%
mutate(
before_flat = list(timeline_points$before_flat),
flat = list(timeline_points$flat),
after_flat = list(timeline_points$after_flat),
delta_t1 = timeline_points$delta_t1,
delta_t2 = timeline_points$delta_t2,
delta_t3 = timeline_points$delta_t3,
TOE = timeline_points$TOE
) %>%
select(-timeline_points) %>%
ungroup()
final_results <- final_results %>%
mutate(
before_flat = sapply(before_flat, function(x) toString(x$Time)),
flat = sapply(flat, function(x) toString(x$Time)),
after_flat = sapply(after_flat, function(x) toString(x$Time))
)
# Save the final results as an Excel file
write_xlsx(final_results, paste0("timeline_results__ethanol_water", temperature, ".xlsx"))
mean_sd_combinations_table <- final_results %>%
group_by(condition) %>%
do(calculate_mean_sd_combinations(.)) %>%
ungroup()
# Find the combination with the smallest combined SD for delta_t1, delta_t2, delta_t3, and TOE
best_combinations_table <- mean_sd_combinations_table %>%
group_by(condition) %>%
mutate(total_sd = sd_delta_t1 + sd_delta_t2 + sd_delta_t3 + sd_toe) %>%
slice_min(total_sd) %>%
ungroup()
# Apply the function to calculate mean Cap for each row of best_combinations_table
best_combinations_table <- best_combinations_table %>%
rowwise() %>%
mutate(
mean_cap = calculate_mean_cap(data, as.character(condition), as.character(combination))$mean_cap,
sd_cap = calculate_mean_cap(data, as.character(condition), as.character(combination))$sd_cap
) %>%
ungroup()
# Define the correct order of conditions
condition_order <- c(
"Pure Water", "10% Ethanol", "20% Ethanol", "30% Ethanol",
"40% Ethanol", "50% Ethanol", "60% Ethanol", "70% Ethanol",
"80% Ethanol", "90% Ethanol", "Pure Ethanol"
)
# Convert the condition column to a factor with the specified levels
best_combinations_table$condition <- factor(best_combinations_table$condition, levels = condition_order)
# Arrange the data frame based on the factor levels
best_combinations_table <- best_combinations_table %>%
arrange(condition)
# Save the tables to Excel
write_xlsx(list(
Mean_SD_Combinations = mean_sd_combinations_table,
Best_Combinations = best_combinations_table
), paste0("mean_sd_combinations_and_selected_runs_ethanol_water", temperature, ".xlsx"))
return(best_combinations_table)
}
# Process for each temperature
data_25C <- read.xlsx("reshape_Water-Ethanol.xlsx", sheet = "25C", colNames = TRUE)
result_25C <- process_data(data_25C, "25C")
data_40C <- read.xlsx("reshape_Water-Ethanol.xlsx", sheet = "40C", colNames = TRUE)
result_40C <- process_data(data_40C, "40C")
data_50C <- read.xlsx("reshape_Water-Ethanol.xlsx", sheet = "50C", colNames = TRUE)
result_50C <- process_data(data_50C, "50C")
data_60C <- read.xlsx("reshape_Water-Ethanol.xlsx", sheet = "60C", colNames = TRUE)
result_60C <- process_data(data_60C, "60C")
##################boxplot###################
combined_results <- bind_rows(
result_25C %>% mutate(temperature = "25C"),
result_40C %>% mutate(temperature = "40C"),
result_50C %>% mutate(temperature = "50C"),
result_60C %>% mutate(temperature = "60C")
)
condition_order <- c(
"Pure Water", "10% Ethanol", "20% Ethanol", "30% Ethanol",
"40% Ethanol", "50% Ethanol", "60% Ethanol", "70% Ethanol",
"80% Ethanol", "90% Ethanol", "Pure Ethanol"
)
# Replace 'Pure Water' with '0' and remove 'Ethanol' and '%' from other labels
clean_condition_labels <- gsub("Pure Water", "0", gsub(" Ethanol|%", "", condition_order))
# Create the plot with summarized data and customizations
combined_results <- combined_results %>%
mutate(Condition_Num = as.numeric(factor(condition, levels = condition_order)))
# Step 1: Summarize the data to calculate the overall mean and SD for each temperature
temperature_summary <- combined_results %>%
group_by(temperature) %>%
summarize(
overall_mean_toe = mean(mean_delta_t2),
overall_sd_toe = sd(mean_delta_t2),
.groups = 'drop'
)
# Step 2: Create the plot with the summarized data
p <- ggplot(combined_results, aes(x = factor(temperature, levels = c("25C", "40C", "50C", "60C")), y = mean_delta_t2)) +
geom_boxplot(outlier.shape = NA) +
geom_errorbar(data = temperature_summary,
aes(x = temperature, ymin = overall_mean_toe - overall_sd_toe, ymax = overall_mean_toe + overall_sd_toe),
width = 0.2, linewidth = 1, color = "black", inherit.aes = FALSE) + # Single error bar for each temperature
geom_jitter(aes(color = factor(Condition_Num)), size = 4, width = 0.1) +
# Use dark colors for each condition
scale_color_manual(
values = c("navyblue", "brown", "#B8860B", "darkgreen", "darkred", "#8B4513", "#4682B4", "#556B2F", "#8B008B", "#2E8B57", "#483D8B"),
breaks = seq_along(condition_order),
labels = clean_condition_labels, # Updated labels with '0' for 'Pure Water' and cleaned others
guide = guide_legend(ncol = 2) # Make legend two columns
) +
# Update axis labels and legend title
labs(
x = "Temperature",
y = "TOE (s)",
color = "Ethanol (%)"
) +
# Customize the theme
theme_minimal() +
theme(
axis.text= element_text(size = 18, family = "Times New Roman",face="bold"),
axis.title = element_text(size = 18, family = "Times New Roman",face="bold"),
legend.text = element_text(size = 18, family = "Times New Roman",face="bold"),
legend.title = element_text(size = 18, family = "Times New Roman",face="bold"),
legend.position = c(0.85, 0.5),
legend.direction = "vertical",
panel.background = element_blank(),
plot.background = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.border = element_rect(color = "black", fill = NA)
)
# Add statistical comparisons using the Wilcoxon test
p <- p + stat_compare_means(comparisons = list(
c("25C", "40C"), c("25C", "50C"),
c("25C", "60C"), c("40C", "50C"),
c("40C", "60C"), c("50C", "60C")
),
method = "t.test", method.args = list(var.equal = FALSE), size = 3.5, family = "Times New Roman")
ggsave("BOXPLOT_t2_E.png", plot = p, width =8, height = 6)
##############RMSE_t2###########
# Combine results for all temperatures
combined_results <- bind_rows(
result_25C %>% mutate(temperature = "25C"),
result_40C %>% mutate(temperature = "40C"),
result_50C %>% mutate(temperature = "50C"),
result_60C %>% mutate(temperature = "60C")
)
# Define condition order for the factor
condition_order <- c(
"Pure Water", "10% Ethanol", "20% Ethanol", "30% Ethanol",
"40% Ethanol", "50% Ethanol", "60% Ethanol", "70% Ethanol",
"80% Ethanol", "90% Ethanol", "Pure Ethanol"
)
combined_results$condition <- factor(combined_results$condition, levels = condition_order)
# Assign to temp_data_t2
temp_data_t2 <- combined_results
# Define custom colors for each temperature
custom_colors <- c("25C" = "navyblue", "40C" = "darkgreen", "50C" = "#B8860B", "60C" = "darkred")
# Add numeric_condition for easier modeling and plotting
temp_data_t2 <- temp_data_t2 %>%
mutate(
numeric_condition = case_when(
condition == "Pure Water" ~ 0,
condition == "Pure Ethanol" ~ 100,
TRUE ~ as.numeric(gsub("[^0-9]", "", as.character(condition)))
)
)
# Ensure a continuous range of numeric_condition
temp_data_t2 <- temp_data_t2 %>%
arrange(temperature, numeric_condition) %>%
group_by(temperature) %>%
complete(numeric_condition = full_seq(numeric_condition, 1)) %>%
ungroup()
# Fit a single linear model after completing missing rows
lm_model <- lm(mean_delta_t2 ~ numeric_condition * as.factor(temperature),
data = temp_data_t2 %>% filter(!is.na(mean_delta_t2)))
# Add linear predicted values
temp_data_t2 <- temp_data_t2 %>%
mutate(predicted_values = predict(lm_model, newdata = temp_data_t2))
# Fit LOESS model per temperature
temp_data_t2 <- temp_data_t2 %>%
group_by(temperature) %>%
mutate(
loess_fit = list(loess(mean_delta_t2 ~ numeric_condition, data = cur_data(), span = 0.75)),
predicted_loess = predict(loess_fit[[1]], newdata = cur_data())
) %>%
ungroup()
# Find min and max values for plotting
min_value_t2 <- min(temp_data_t2$mean_delta_t2, na.rm = TRUE)
max_value_t2 <- max(temp_data_t2$mean_delta_t2, na.rm = TRUE)
# Calculate accuracy metrics for linear model (using a cubic polynomial model for evaluation)
temp_models_t2 <- temp_data_t2 %>%
filter(!is.na(mean_delta_t2)) %>%
group_by(temperature) %>%
nest() %>%
mutate(
model = map(data, ~ lm(mean_delta_t2 ~ poly(numeric_condition, 3), data = .x)),
r_squared = map_dbl(model, ~ glance(.x)$r.squared),
rmse = map_dbl(model, ~ sqrt(mean(residuals(.x)^2))),
mae = map_dbl(model, ~ mean(abs(residuals(.x)))),
rmse_percentage = map_dbl(data, ~ {
preds <- predict(lm(mean_delta_t2 ~ poly(numeric_condition, 3), data = .x))
range_y <- max(.x$mean_delta_t2) - min(.x$mean_delta_t2)
rmse_val <- sqrt(mean((.x$mean_delta_t2 - preds)^2))
(rmse_val / range_y) * 100
})
)
# Calculate accuracy metrics for LOESS model
temp_models_t2_loess <- temp_data_t2 %>%
filter(!is.na(mean_delta_t2)) %>%
group_by(temperature) %>%
summarize(
r_squared_loess = {
loess_resid = mean_delta_t2 - predicted_loess
1 - (sum(loess_resid^2) / sum((mean_delta_t2 - mean(mean_delta_t2))^2))
},
rmse_loess = {
loess_resid = mean_delta_t2 - predicted_loess
sqrt(mean(loess_resid^2))
},
mae_loess = {
loess_resid = mean_delta_t2 - predicted_loess
mean(abs(loess_resid))
},
rmse_percentage_loess = {
loess_resid = mean_delta_t2 - predicted_loess
range_y <- max(mean_delta_t2) - min(mean_delta_t2)
(sqrt(mean(loess_resid^2)) / range_y) * 100
}
)
# Combine accuracy metrics
accuracy_combined <- temp_models_t2 %>%
select(temperature, rmse_percentage) %>%
left_join(
temp_models_t2_loess %>% select(temperature, rmse_percentage_loess),
by = "temperature"
)
# Combine accuracy metrics
accuracy_combined <- temp_models_t2 %>%
select(temperature, rmse_percentage) %>%
left_join(
temp_models_t2_loess %>% select(temperature, rmse_percentage_loess),
by = "temperature"
)
# Prepare a table with 4 rows and 2 columns for annotation
my_table <- data.frame(
Temperature = accuracy_combined$temperature,
`RMSE (%) (Linear)` = round(accuracy_combined$rmse_percentage, 3),
`RMSE (%) (LOESS)` = round(accuracy_combined$rmse_percentage_loess, 3)
)
colnames(my_table) <- c("Temperature", "RMSE_Linear(%)", "RMSE_LOESS(%)")
tt <- ttheme_minimal(
core = list(
fg_params = list(
fontface = c("bold", "bold", "bold", "bold"), # Bold text for all rows
col = c("navyblue", "darkgreen", "#B8860B", "darkred") # Row colors
),
bg_params = list(col = "black", lwd = 1)
),
colhead = list(
fg_params = list(fontface = "bold", col = "black"), # Black title
bg_params = list(col = "black", lwd = 1)
)
)
# Apply the custom theme to your table
table_grob <- tableGrob(my_table, rows = NULL, theme = tt)
# Split data into segments for numeric_condition <= 80 and > 80
temp_data_t2$segment <- ifelse(temp_data_t2$numeric_condition > 80, "above_80", "below_80")
p_t2 <- ggplot(temp_data_t2, aes(x = numeric_condition, y = mean_delta_t2, color = as.factor(temperature))) +
geom_point(size = 3, na.rm = TRUE) + # Points for original data
# LOESS dashed line
geom_line(aes(y = predicted_loess, group = temperature, linetype = "LOESS"), size = 1) +
# Linear solid line for numeric_condition <= 80
geom_line(
data = temp_data_t2 %>% filter(segment == "below_80"),
aes(y = predicted_values, group = temperature, linetype = "Linear (under 80%)"),
size = 1.8
) +
# Linear dotted line for numeric_condition > 80
geom_line(
data = temp_data_t2 %>% filter(segment == "above_80"),
aes(y = predicted_values, group = temperature, linetype = "Linear (80%-100%)"),
size = 1.2
) +
labs(
x = "Ethanol (%)",
y = expression(bold("Mean"~Delta~t[2]~"(s)")),
linetype = "Line Type"
) +
theme_minimal() +
theme(
axis.text.x = element_text(size = 18, family = "Times New Roman", face = "bold"),
axis.text.y = element_text(size = 18, family = "Times New Roman", face = "bold"),
axis.title = element_text(size = 18, family = "Times New Roman", face = "bold"),
legend.text = element_text(size = 15, family = "Times New Roman", face = "bold"),
legend.title = element_text(size = 15, family = "Times New Roman", face = "bold"),
legend.key.height = unit(1.5, "lines"),
legend.position = c(0.83, 0.55),
panel.background = element_blank(),
plot.background = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.border = element_rect(color = "black", fill = NA)
) +
scale_color_manual(values = custom_colors, guide = "none") + # Remove color legend
scale_linetype_manual(
values = c(
"Linear (under 80%)" = "solid",
"Linear (80%-100%)" = "dotted",
"LOESS" = "dashed"
),
breaks = c("Linear (under 80%)", "Linear (80%-100%)", "LOESS") # Set the order explicitly
) +
coord_cartesian(ylim = c(min_value_t2, max_value_t2)) +
scale_x_continuous(breaks = seq(0, 100, 10), labels = as.character(seq(0, 100, 10))) +
# Add the table annotation
annotation_custom(
grob = table_grob,
xmin = 50, xmax = 90, # Adjust as needed for positioning
ymin = max_value_t2 - 100,
ymax = max_value_t2
)
# Save the plot
ggsave("RMSE_t2_E.png", plot = p_t2, width = 8, height = 6)
##############RMSE_cap###########
# Define condition order for proper factor levels
condition_order <- c(
"Pure Water", "10% Ethanol", "20% Ethanol", "30% Ethanol",
"40% Ethanol", "50% Ethanol", "60% Ethanol", "70% Ethanol",
"80% Ethanol", "90% Ethanol", "Pure Ethanol"
)
# Combine results for all temperatures
all_results <- bind_rows(
result_25C %>% mutate(Temperature = "25C"),
result_40C %>% mutate(Temperature = "40C"),
result_50C %>% mutate(Temperature = "50C"),
result_60C %>% mutate(Temperature = "60C")
)
# Convert 'condition' to factor with specified levels
all_results$condition <- factor(all_results$condition, levels = condition_order)
# Summarize data by condition and temperature
all_results_summarized <- all_results %>%
group_by(condition, Temperature) %>%
summarize(
mean_cap = mean(mean_cap),
sd_cap = mean(sd_cap),
.groups = 'drop'
)
# Add numeric_condition for easier modeling
all_results_summarized <- all_results_summarized %>%
mutate(numeric_condition = as.numeric(gsub("[^0-9]", "", as.character(condition))),
numeric_condition = ifelse(is.na(numeric_condition),
ifelse(condition == "Pure Water", 0, 100), numeric_condition))
# Fit models for mean_cap using both Linear and LOESS
temp_models_cap <- all_results_summarized %>%
group_by(Temperature) %>%
nest() %>%
mutate(
# Fit Linear model
linear_model = map(data, ~ lm(mean_cap ~ poly(numeric_condition, 3), data = .x)),
# Fit LOESS model
loess_model = map(data, ~ loess(mean_cap ~ numeric_condition, data = .x, span = 0.75)),
# Calculate predictions for Linear and LOESS
linear_predictions = map(linear_model, ~ augment(.x)),
loess_predictions = map(loess_model, ~ data.frame(.x$fitted)),
# Calculate accuracy metrics
r_squared_linear = map_dbl(linear_model, ~ glance(.x)$r.squared),
rmse_linear = map_dbl(linear_model, ~ sqrt(mean(residuals(.x)^2))),
rmse_loess = map_dbl(loess_model, ~ sqrt(mean((.x$residuals)^2))),
mae_linear = map_dbl(linear_model, ~ mean(abs(residuals(.x))))
) %>%
mutate(
range_y = map_dbl(data, ~ max(.x$mean_cap) - min(.x$mean_cap)), # Range of the observed values
rmse_percentage_linear = (rmse_linear / range_y) * 100,
rmse_percentage_loess = (rmse_loess / range_y) * 100
)
# Print accuracy metrics for each temperature
temp_models_cap %>%
select(Temperature, r_squared_linear, rmse_percentage_linear, rmse_percentage_loess)
## # A tibble: 4 × 4
## # Groups: Temperature [4]
## Temperature r_squared_linear rmse_percentage_linear rmse_percentage_loess
## <chr> <dbl> <dbl> <dbl>
## 1 25C 0.903 9.56 6.30
## 2 40C 0.826 12.9 12.9
## 3 50C 0.885 10.1 7.05
## 4 60C 0.781 14.3 10.8
# Find min and max values for the plot
min_value_cap <- min(all_results_summarized$mean_cap, na.rm = TRUE)
max_value_cap <- max(all_results_summarized$mean_cap, na.rm = TRUE) + 5 # Add buffer to max for annotations
# Prepare a table with RMSE for both Linear and LOESS
accuracy_combined <- temp_models_cap %>%
select(
Temperature,
`RMSE% (Linear)` = rmse_percentage_linear,
`RMSE% (LOESS)` = rmse_percentage_loess
)
# Define custom theme for table
tt <- ttheme_minimal(
core = list(
fg_params = list(
fontface = c("bold", "bold", "bold", "bold"), # Bold text for all rows
col = c("navyblue", "darkgreen", "#B8860B", "darkred") # Row colors
),
bg_params = list(col = "black", lwd = 1)
),
colhead = list(
fg_params = list(fontface = "bold", col = "black"), # Black title
bg_params = list(col = "black", lwd = 1)
)
)
# Create table grob
table_grob <- tableGrob(accuracy_combined, rows = NULL, theme = tt)
# Add a separate legend for line types
p <- ggplot(all_results_summarized, aes(x = numeric_condition, y = mean_cap, group = Temperature)) +
geom_point(aes(color = Temperature), size = 3, show.legend = FALSE) + # Points for data, no legend
geom_smooth(aes(linetype = "Linear", color = Temperature), method = "lm", se = FALSE, size = 1) + # Linear model
geom_smooth(aes(linetype = "LOESS", color = Temperature), method = "loess", se = FALSE, size = 1.8) + # LOESS model
labs(x = "Ethanol (%)", y = expression(bold(Mean~Delta~Cap~(fF))), linetype = "Model") +
theme_minimal() +
theme(
axis.text.x = element_text(size = 18, family = "Times New Roman", face = "bold"),
axis.text.y = element_text(size = 18, family = "Times New Roman", face = "bold"),
axis.title = element_text(size = 18, family = "Times New Roman", face = "bold"),
legend.text = element_text(size = 18, family = "Times New Roman", face = "bold"),
legend.title = element_text(size = 18, family = "Times New Roman", face = "bold"),
panel.background = element_blank(),
plot.background = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.border = element_rect(color = "black", fill = NA),
legend.position = c(0.1, 0.15) # Adjust legend position as needed
) +
coord_cartesian(ylim = c(min_value_cap, max_value_cap)) +
scale_color_manual(values = c("25C" = "navyblue", "40C" = "darkgreen", "50C" = "#B8860B", "60C" = "darkred"), guide = "none") +
scale_linetype_manual(
values = c("Linear" = "solid", "LOESS" = "dashed"),
breaks = c("Linear", "LOESS") # Order of the legend
) +
guides(
linetype = guide_legend(override.aes = list(color = "black")) # Black line symbols in legend
) +
# Add table annotation
annotation_custom(
grob = table_grob,
xmin = 55, xmax = 90, # Adjust as needed for positioning
ymin = max_value_cap - 7,
ymax = max_value_cap
)
# Save the plot
ggsave("RMSE_cap_E.png", plot = p, width = 8, height = 6)
#_______________________Metanol_______________________-
process_data <- function(data, temperature) {
# Replace "-run" with "_run" in the Condition column to ensure uniformity
data <- data %>%
mutate(Condition = str_replace(Condition, "-run", "_run"))
unique_conditions <- unique(data$Condition)
split_conditions <- data.frame(
condition = sapply(strsplit(unique_conditions, "_"), `[`, 1),
run_name = sapply(strsplit(unique_conditions, "_"), `[`, 2)
)
results <- split_conditions %>%
mutate(timeline_points = map2(condition, run_name, ~get_timeline_points(data, .x, .y)))
final_results <- results %>%
rowwise() %>%
mutate(
before_flat = list(timeline_points$before_flat),
flat = list(timeline_points$flat),
after_flat = list(timeline_points$after_flat),
delta_t1 = timeline_points$delta_t1,
delta_t2 = timeline_points$delta_t2,
delta_t3 = timeline_points$delta_t3,
TOE = timeline_points$TOE
) %>%
select(-timeline_points) %>%
ungroup()
final_results <- final_results %>%
mutate(
before_flat = sapply(before_flat, function(x) toString(x$Time)),
flat = sapply(flat, function(x) toString(x$Time)),
after_flat = sapply(after_flat, function(x) toString(x$Time))
)
# Save the final results as an Excel file
write_xlsx(final_results, paste0("timeline_results__Methanol_water", temperature, ".xlsx"))
mean_sd_combinations_table <- final_results %>%
group_by(condition) %>%
do(calculate_mean_sd_combinations(.)) %>%
ungroup()
# Find the combination with the smallest combined SD for delta_t1, delta_t2, delta_t3, and TOE
best_combinations_table <- mean_sd_combinations_table %>%
group_by(condition) %>%
mutate(total_sd = sd_delta_t1 + sd_delta_t2 + sd_delta_t3 + sd_toe) %>%
slice_min(total_sd) %>%
ungroup()
# Apply the function to calculate mean Cap for each row of best_combinations_table
best_combinations_table <- best_combinations_table %>%
rowwise() %>%
mutate(
mean_cap = calculate_mean_cap(data, as.character(condition), as.character(combination))$mean_cap,
sd_cap = calculate_mean_cap(data, as.character(condition), as.character(combination))$sd_cap
) %>%
ungroup()
# Define the correct order of conditions
condition_order <- c(
"Pure Water", "10% Methanol", "20% Methanol", "30% Methanol",
"40% Methanol", "50% Methanol", "60% Methanol", "70% Methanol",
"80% Methanol", "90% Methanol", "Pure Methanol"
)
# Convert the condition column to a factor with the specified levels
best_combinations_table$condition <- factor(best_combinations_table$condition, levels = condition_order)
# Arrange the data frame based on the factor levels
best_combinations_table <- best_combinations_table %>%
arrange(condition)
# Save the tables to Excel
write_xlsx(list(
Mean_SD_Combinations = mean_sd_combinations_table,
Best_Combinations = best_combinations_table
), paste0("mean_sd_combinations_and_selected_runs_Methanol_water", temperature, ".xlsx"))
return(best_combinations_table)
}
# Process for each temperature
data_25C <- read.xlsx("reshape_Water-Methanol.xlsx", sheet = "25C", colNames = TRUE)
result_25C <- process_data(data_25C, "25C")
data_40C <- read.xlsx("reshape_Water-Methanol.xlsx", sheet = "40C", colNames = TRUE)
result_40C <- process_data(data_40C, "40C")
data_50C <- read.xlsx("reshape_Water-Methanol.xlsx", sheet = "50C", colNames = TRUE)
result_50C <- process_data(data_50C, "50C")
data_60C <- read.xlsx("reshape_Water-Methanol.xlsx", sheet = "60C", colNames = TRUE)
result_60C <- process_data(data_60C, "60C")
#############boxplot##############
combined_results <- bind_rows(
result_25C %>% mutate(temperature = "25C"),
result_40C %>% mutate(temperature = "40C"),
result_50C %>% mutate(temperature = "50C"),
result_60C %>% mutate(temperature = "60C")
)
condition_order <- c(
"Pure Water", "10% Methanol", "20% Methanol", "30% Methanol",
"40% Methanol", "50% Methanol", "60% Methanol", "70% Methanol",
"80% Methanol", "90% Methanol", "Pure Methanol"
)
# Replace 'Pure Water' with '0' and remove 'Methanol' and '%' from other labels
clean_condition_labels <- gsub("Pure Water", "0", gsub(" Methanol|%", "", condition_order))
# Create the plot with summarized data and customizations
combined_results <- combined_results %>%
mutate(Condition_Num = as.numeric(factor(condition, levels = condition_order)))
# Step 1: Summarize the data to calculate the overall mean and SD for each temperature
temperature_summary <- combined_results %>%
group_by(temperature) %>%
summarize(
overall_mean_toe = mean(mean_toe),
overall_sd_toe = sd(mean_toe),
.groups = 'drop'
)
# Step 2: Create the plot with the summarized data
p <- ggplot(combined_results, aes(x = factor(temperature, levels = c("25C", "40C", "50C", "60C")), y = mean_toe)) +
geom_boxplot(outlier.shape = NA) +
geom_errorbar(data = temperature_summary,
aes(x = temperature, ymin = overall_mean_toe - overall_sd_toe, ymax = overall_mean_toe + overall_sd_toe),
width = 0.2, linewidth = 1, color = "black", inherit.aes = FALSE) + # Single error bar for each temperature
geom_jitter(aes(color = factor(Condition_Num)), size = 4, width = 0.1) +
# Use dark colors for each condition
scale_color_manual(
values = c("navyblue", "brown", "#B8860B", "darkgreen", "darkred", "#8B4513", "#4682B4", "#556B2F", "#8B008B", "#2E8B57", "#483D8B"),
breaks = seq_along(condition_order),
labels = clean_condition_labels, # Updated labels with '0' for 'Pure Water' and cleaned others
guide = guide_legend(ncol = 2) # Make legend two columns
) +
# Update axis labels and legend title
labs(
x = "Temperature",
y = "TOE (s)",
color = "Methanol (%)"
) +
# Customize the theme
theme_minimal() +
theme(
axis.text= element_text(size = 18, family = "Times New Roman",face="bold"),
axis.title = element_text(size = 18, family = "Times New Roman",face="bold"),
legend.text = element_text(size = 18, family = "Times New Roman",face="bold"),
legend.title = element_text(size = 18, family = "Times New Roman",face="bold"),
legend.position = c(0.85, 0.5),
legend.direction = "vertical",
panel.background = element_blank(),
plot.background = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.border = element_rect(color = "black", fill = NA)
)
# Add statistical comparisons using the Wilcoxon test
p <- p + stat_compare_means(comparisons = list(
c("25C", "40C"), c("25C", "50C"),
c("25C", "60C"), c("40C", "50C"),
c("40C", "60C"), c("50C", "60C")
),
method = "t.test", method.args = list(var.equal = FALSE), size = 3.5, family = "Times New Roman",fontface = "bold")
ggsave("boxplot_M.png", plot = p, width =8, height = 6)
#############RMSEt2#################
# Combine results for all temperatures
combined_results <- bind_rows(
result_25C %>% mutate(temperature = "25C"),
result_40C %>% mutate(temperature = "40C"),
result_50C %>% mutate(temperature = "50C"),
result_60C %>% mutate(temperature = "60C")
)
# Define condition order for the factor
condition_order <- c(
"Pure Water", "10% Methanol", "20% Methanol", "30% Methanol",
"40% Methanol", "50% Methanol", "60% Methanol", "70% Methanol",
"80% Methanol", "90% Methanol", "Pure Methanol"
)
combined_results$condition <- factor(combined_results$condition, levels = condition_order)
# Assign to temp_data_t2
temp_data_t2 <- combined_results
# Define custom colors for each temperature
custom_colors <- c("25C" = "navyblue", "40C" = "darkgreen", "50C" = "#B8860B", "60C" = "darkred")
# Add numeric_condition for easier modeling and plotting
temp_data_t2 <- temp_data_t2 %>%
mutate(
numeric_condition = case_when(
condition == "Pure Water" ~ 0,
condition == "Pure Methanol" ~ 100,
TRUE ~ as.numeric(gsub("[^0-9]", "", as.character(condition)))
)
)
# Ensure a continuous range of numeric_condition
temp_data_t2 <- temp_data_t2 %>%
arrange(temperature, numeric_condition) %>%
group_by(temperature) %>%
complete(numeric_condition = full_seq(numeric_condition, 1)) %>%
ungroup()
# Fit a single linear model after completing missing rows
lm_model <- lm(mean_delta_t2 ~ numeric_condition * as.factor(temperature),
data = temp_data_t2 %>% filter(!is.na(mean_delta_t2)))
# Add linear predicted values
temp_data_t2 <- temp_data_t2 %>%
mutate(predicted_values = predict(lm_model, newdata = temp_data_t2))
# Fit LOESS model per temperature
temp_data_t2 <- temp_data_t2 %>%
group_by(temperature) %>%
mutate(
loess_fit = list(loess(mean_delta_t2 ~ numeric_condition, data = cur_data(), span = 0.75)),
predicted_loess = predict(loess_fit[[1]], newdata = cur_data())
) %>%
ungroup()
# Find min and max values for plotting
min_value_t2 <- min(temp_data_t2$mean_delta_t2, na.rm = TRUE)
max_value_t2 <- max(temp_data_t2$mean_delta_t2, na.rm = TRUE)
# Calculate accuracy metrics for linear model (using a cubic polynomial model for evaluation)
temp_models_t2 <- temp_data_t2 %>%
filter(!is.na(mean_delta_t2)) %>%
group_by(temperature) %>%
nest() %>%
mutate(
model = map(data, ~ lm(mean_delta_t2 ~ poly(numeric_condition, 3), data = .x)),
r_squared = map_dbl(model, ~ glance(.x)$r.squared),
rmse = map_dbl(model, ~ sqrt(mean(residuals(.x)^2))),
mae = map_dbl(model, ~ mean(abs(residuals(.x)))),
rmse_percentage = map_dbl(data, ~ {
preds <- predict(lm(mean_delta_t2 ~ poly(numeric_condition, 3), data = .x))
range_y <- max(.x$mean_delta_t2) - min(.x$mean_delta_t2)
rmse_val <- sqrt(mean((.x$mean_delta_t2 - preds)^2))
(rmse_val / range_y) * 100
})
)
# Calculate accuracy metrics for LOESS model
temp_models_t2_loess <- temp_data_t2 %>%
filter(!is.na(mean_delta_t2)) %>%
group_by(temperature) %>%
summarize(
r_squared_loess = {
loess_resid = mean_delta_t2 - predicted_loess
1 - (sum(loess_resid^2) / sum((mean_delta_t2 - mean(mean_delta_t2))^2))
},
rmse_loess = {
loess_resid = mean_delta_t2 - predicted_loess
sqrt(mean(loess_resid^2))
},
mae_loess = {
loess_resid = mean_delta_t2 - predicted_loess
mean(abs(loess_resid))
},
rmse_percentage_loess = {
loess_resid = mean_delta_t2 - predicted_loess
range_y <- max(mean_delta_t2) - min(mean_delta_t2)
(sqrt(mean(loess_resid^2)) / range_y) * 100
}
)
# Combine accuracy metrics
accuracy_combined <- temp_models_t2 %>%
select(temperature, rmse_percentage) %>%
left_join(
temp_models_t2_loess %>% select(temperature, rmse_percentage_loess),
by = "temperature"
)
# Combine accuracy metrics
accuracy_combined <- temp_models_t2 %>%
select(temperature, rmse_percentage) %>%
left_join(
temp_models_t2_loess %>% select(temperature, rmse_percentage_loess),
by = "temperature"
)
# Prepare a table with 4 rows and 2 columns for annotation
my_table <- data.frame(
Temperature = accuracy_combined$temperature,
`RMSE (%) (Linear)` = round(accuracy_combined$rmse_percentage, 3),
`RMSE (%) (LOESS)` = round(accuracy_combined$rmse_percentage_loess, 3)
)
colnames(my_table) <- c("Temperature", "RMSE_Linear(%)", "RMSE_LOESS(%)")
tt <- ttheme_minimal(
core = list(
fg_params = list(
fontface = c("bold", "bold", "bold", "bold"), # Bold text for all rows
col = c("navyblue", "darkgreen", "#B8860B", "darkred") # Row colors
),
bg_params = list(col = "black", lwd = 1)
),
colhead = list(
fg_params = list(fontface = "bold", col = "black"), # Black title
bg_params = list(col = "black", lwd = 1)
)
)
# Apply the custom theme to your table
table_grob <- tableGrob(my_table, rows = NULL, theme = tt)
# Split data into segments for numeric_condition <= 80 and > 80
temp_data_t2$segment <- ifelse(temp_data_t2$numeric_condition > 80, "above_80", "below_80")
p_t2 <- ggplot(temp_data_t2, aes(x = numeric_condition, y = mean_delta_t2, color = as.factor(temperature))) +
geom_point(size = 3, na.rm = TRUE) + # Points for original data
# LOESS dashed line
geom_line(aes(y = predicted_loess, group = temperature, linetype = "LOESS"), size = 1) +
# Linear solid line for numeric_condition <= 80
geom_line(
data = temp_data_t2 %>% filter(segment == "below_80"),
aes(y = predicted_values, group = temperature, linetype = "Linear (under 80%)"),
size = 1.8
) +
# Linear dotted line for numeric_condition > 80
geom_line(
data = temp_data_t2 %>% filter(segment == "above_80"),
aes(y = predicted_values, group = temperature, linetype = "Linear (80%-100%)"),
size = 1.2
) +
labs(
x = "Methanol (%)",
y = expression(bold("Mean"~Delta~t[2]~"(s)")),
linetype = "Line Type"
) +
theme_minimal() +
theme(
axis.text.x = element_text(size = 18, family = "Times New Roman", face = "bold"),
axis.text.y = element_text(size = 18, family = "Times New Roman", face = "bold"),
axis.title = element_text(size = 18, family = "Times New Roman", face = "bold"),
legend.text = element_text(size = 15, family = "Times New Roman", face = "bold"),
legend.title = element_text(size = 15, family = "Times New Roman", face = "bold"),
legend.key.height = unit(1.5, "lines"),
legend.position = c(0.83, 0.55),
panel.background = element_blank(),
plot.background = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.border = element_rect(color = "black", fill = NA)
) +
scale_color_manual(values = custom_colors, guide = "none") + # Remove color legend
scale_linetype_manual(
values = c(
"Linear (under 80%)" = "solid",
"Linear (80%-100%)" = "dotted",
"LOESS" = "dashed"
),
breaks = c("Linear (under 80%)", "Linear (80%-100%)", "LOESS") # Set the order explicitly
) +
coord_cartesian(ylim = c(min_value_t2, max_value_t2)) +
scale_x_continuous(breaks = seq(0, 100, 10), labels = as.character(seq(0, 100, 10))) +
# Add the table annotation
annotation_custom(
grob = table_grob,
xmin = 50, xmax = 90, # Adjust as needed for positioning
ymin = max_value_t2 - 100,
ymax = max_value_t2
)
ggsave("RMSE_t2_M.png", plot = p_t2, width =8, height = 6)
############rmse_cap#############
###################
condition_order <- c(
"Pure Water", "10% Methanol", "20% Methanol", "30% Methanol",
"40% Methanol", "50% Methanol", "60% Methanol", "70% Methanol",
"80% Methanol", "90% Methanol", "Pure Methanol"
)
# Combine results for all temperatures
all_results <- bind_rows(
result_25C %>% mutate(Temperature = "25C"),
result_40C %>% mutate(Temperature = "40C"),
result_50C %>% mutate(Temperature = "50C"),
result_60C %>% mutate(Temperature = "60C")
)
# Convert 'condition' to factor with specified levels
all_results$condition <- factor(all_results$condition, levels = condition_order)
# Summarize data by condition and temperature
all_results_summarized <- all_results %>%
group_by(condition, Temperature) %>%
summarize(
mean_cap = mean(mean_cap),
sd_cap = mean(sd_cap),
.groups = 'drop'
)
# Add numeric_condition for easier modeling
all_results_summarized <- all_results_summarized %>%
mutate(numeric_condition = as.numeric(gsub("[^0-9]", "", as.character(condition))),
numeric_condition = ifelse(is.na(numeric_condition),
ifelse(condition == "Pure Water", 0, 100), numeric_condition))
# Fit models for mean_cap using both Linear and LOESS
temp_models_cap <- all_results_summarized %>%
group_by(Temperature) %>%
nest() %>%
mutate(
# Fit Linear model
linear_model = map(data, ~ lm(mean_cap ~ poly(numeric_condition, 3), data = .x)),
# Fit LOESS model
loess_model = map(data, ~ loess(mean_cap ~ numeric_condition, data = .x, span = 0.75)),
# Calculate predictions for Linear and LOESS
linear_predictions = map(linear_model, ~ augment(.x)),
loess_predictions = map(loess_model, ~ data.frame(.x$fitted)),
# Calculate accuracy metrics
r_squared_linear = map_dbl(linear_model, ~ glance(.x)$r.squared),
rmse_linear = map_dbl(linear_model, ~ sqrt(mean(residuals(.x)^2))),
rmse_loess = map_dbl(loess_model, ~ sqrt(mean((.x$residuals)^2))),
mae_linear = map_dbl(linear_model, ~ mean(abs(residuals(.x))))
) %>%
mutate(
range_y = map_dbl(data, ~ max(.x$mean_cap) - min(.x$mean_cap)), # Range of the observed values
rmse_percentage_linear = (rmse_linear / range_y) * 100,
rmse_percentage_loess = (rmse_loess / range_y) * 100
)
# Print accuracy metrics for each temperature
temp_models_cap %>%
select(Temperature, r_squared_linear, rmse_percentage_linear, rmse_percentage_loess)
## # A tibble: 4 × 4
## # Groups: Temperature [4]
## Temperature r_squared_linear rmse_percentage_linear rmse_percentage_loess
## <chr> <dbl> <dbl> <dbl>
## 1 25C 0.529 18.4 17.2
## 2 40C 0.805 11.1 11.1
## 3 50C 0.794 14.1 12.0
## 4 60C 0.585 22.1 20.4
# Find min and max values for the plot
min_value_cap <- min(all_results_summarized$mean_cap, na.rm = TRUE)
max_value_cap <- max(all_results_summarized$mean_cap, na.rm = TRUE) + 5 # Add buffer to max for annotations
# Prepare a table with RMSE for both Linear and LOESS
accuracy_combined <- temp_models_cap %>%
select(
Temperature,
`RMSE% (Linear)` = rmse_percentage_linear,
`RMSE% (LOESS)` = rmse_percentage_loess
)
# Define custom theme for table
tt <- ttheme_minimal(
core = list(
fg_params = list(
fontface = c("bold", "bold", "bold", "bold"), # Bold text for all rows
col = c("navyblue", "darkgreen", "#B8860B", "darkred") # Row colors
),
bg_params = list(col = "black", lwd = 1)
),
colhead = list(
fg_params = list(fontface = "bold", col = "black"), # Black title
bg_params = list(col = "black", lwd = 1)
)
)
# Create table grob
table_grob <- tableGrob(accuracy_combined, rows = NULL, theme = tt)
# Add a separate legend for line types
p <- ggplot(all_results_summarized, aes(x = numeric_condition, y = mean_cap, group = Temperature)) +
geom_point(aes(color = Temperature), size = 3, show.legend = FALSE) + # Points for data, no legend
geom_smooth(aes(linetype = "Linear", color = Temperature), method = "lm", se = FALSE, size = 1) + # Linear model
geom_smooth(aes(linetype = "LOESS", color = Temperature), method = "loess", se = FALSE, size = 1.8) + # LOESS model
labs(x = "Methanol (%)", y = expression(bold(Mean~Delta~Cap~(fF))), linetype = "Model") +
theme_minimal() +
theme(
axis.text.x = element_text(size = 18, family = "Times New Roman", face = "bold"),
axis.text.y = element_text(size = 18, family = "Times New Roman", face = "bold"),
axis.title = element_text(size = 18, family = "Times New Roman", face = "bold"),
legend.text = element_text(size = 18, family = "Times New Roman", face = "bold"),
legend.title = element_text(size = 18, family = "Times New Roman", face = "bold"),
panel.background = element_blank(),
plot.background = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.border = element_rect(color = "black", fill = NA),
legend.position = c(0.1, 0.15) # Adjust legend position as needed
) +
coord_cartesian(ylim = c(min_value_cap, max_value_cap)) +
scale_color_manual(values = c("25C" = "navyblue", "40C" = "darkgreen", "50C" = "#B8860B", "60C" = "darkred"), guide = "none") +
scale_linetype_manual(
values = c("Linear" = "solid", "LOESS" = "dashed"),
breaks = c("Linear", "LOESS") # Order of the legend
) +
guides(
linetype = guide_legend(override.aes = list(color = "black")) # Black line symbols in legend
) +
# Add table annotation
annotation_custom(
grob = table_grob,
xmin = 55, xmax = 90, # Adjust as needed for positioning
ymin = max_value_cap - 7,
ymax = max_value_cap
)
ggsave("RMSE_cap_M.png", plot = p, width =8, height = 6)
#_______________________-arrange
#___________________________________
library(magick)
# Load and resize the images
RMSE_t2_E <- image_read("RMSE_t2_E.png") %>% image_background("white") %>% image_resize("2000x600") # Make RMSE_t2 bigger
RMSE_cap_E <- image_read("RMSE_cap_E.png") %>% image_background("white") %>% image_resize("800x600") # Shrink RMSE_cap
RMSE_t2_M <- image_read("RMSE_t2_M.png") %>% image_background("white") %>% image_resize("2000x600") # Make RMSE_t2 bigger
RMSE_cap_M <- image_read("RMSE_cap_M.png") %>% image_background("white") %>% image_resize("800x600") # Shrink RMSE_cap to
boxplot_Ethanol <- image_read("BOXPLOT_t2_E.png") %>% image_background("white") %>% image_resize("2000x600") # Make RMSE_t2 bigger
boxplot_Methanol<- image_read("boxplot_M.png") %>% image_background("white") %>% image_resize("800x600") # Shrink RMSE_cap to
# Arrange the first row: RMSE_cap_E, RMSE_t2_E
# Create horizontal spacer
horizontal_spacer <- image_blank(width = 50, height = image_info(RMSE_cap_E)$height, color = "white") # 50px wide blank space
row1 <- image_append(c(RMSE_t2_E, horizontal_spacer, RMSE_t2_M))
row2 <- image_append(c(RMSE_cap_E, horizontal_spacer, RMSE_cap_M))
row3 <- image_append(c(boxplot_Ethanol, horizontal_spacer, boxplot_Methanol))
# Create vertical spacer for alignment
vertical_spacer <- image_blank(width = image_info(row1)$width, height = 50, color = "white") # 200px height blank space
# Stack the rows vertically with spacers between them
final_image <- image_append(c(row1, vertical_spacer, row2, vertical_spacer, row3), stack = TRUE)
# Save the final combined image
image_write(final_image, "combined_plot1.png")
knitr::include_graphics("combined_plot1.png")
